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This  work  quantifies  relationships  between  the  spatial,  or  Eulerian,  distribution  of 
the  properties  of  a  chemically  and  physically  heterogeneous  porous  medium  and  those 
as  observed  along  the  natural,  or  Lagrangian,  trajectories  that  a  fluid  particle  traces 
in  a  steady  and  irrotational  flow.  From  these  relationships,  expressions  that  relate 
the  transport  of  solutes  through  the  porous  medium  along  the  natural  trajectories 
to  the  aforementioned  Eulerian  distributions  are  developed.  The  effects  of  injection 
mode  upon  global  measures  of  transport  as  reflected  by  the  temporal  moments  of 
breakthrough  curves  and  spatial  moments  of  a  solute  plume  are  developed.  The 
coupled  effects  of  correlation  of  a  linear  equilibrium  sorption  to  the  underlying  log 
hydraulic  conductivity  field  and  injection  mode  on  the  evolving  temporal  moments 
of  mass  breakthrough  curve  and  the  coupled  effects  of  correlation  of  a  first-order 
decay  coefficient  and  injection  mode  upon  the  spatial  moments  of  a  solute  plume  are 
examined. 
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CHAPTER  1 
INTRODUCTION 

1.1    Primary  Contribution 

The  primary  contribution  of  this  work  is  documenting  a  consistent  approach  to 
the  development  of  estimators  of  solute  transport  in  the  Lagrangian  framework.  This 
consistent  approach  is  to  identify  the  collection  of  Lagrangian  trajectories  associated 
with  a  stationary  Eulerian  velocity  field  for  which  the  Lagrangian  velocities  are  also 
stationary.  This  is  a  modest  extension  of  work  pioneered  by  myriad  hydrologists. 
The  benefit  derived  from  this  contribution  is  enhanced  understanding  of  transport  of 
solutes  in  heterogeneous  velocity  fields. 

1.2    Historical  Context 

A  grant  to  study  the  impact  of  heterogeneous  source  morphology  upon  subse- 
quent solute  transport  initiated  this  research.  In  seeking  a  suitable  framework  and 
"toolset"  with  which  to  attack  this  problem,  there  appeared  what  seemed  to  be  funda- 
mental inconsistencies  in  the  literature.  Given  the  peculiar  tendencies  of  the  author, 
the  "larger  issues"  could  not  be  addressed  until  these  apparent  inconsistencies  were 
resolved. 

Gedeon  Dagan  developed  a  robust  theory  of  solute  transport  using  a  stochastic- 
Lagrangian  framework  [Dagan,  1982a,  b].  In  fact,  an  entire  "school"  of  stochastic 
subsurface  hydrology  sprung  from  his  work. 

Allen  Shapiro  and  Vladimir  Cvetkovic  wrote  Stochastic  analysis  of  solute  arrival 
time  in  heterogeneous  porous  media  in  1988  [Shapiro  and  Cvetkovic,  1988].  The  au- 
thors made  an  explicit  assumption  often  implicitly  made  in  the  Lagrangian  transport 
literature,  namely  that  a  fluid  parcel  deviates  little  from  its  mean  trajectory  in  weakly 
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heterogeneous  conductivity  fields,  and  that  the  Lagrangian  and  Eulerian  velocities  are 
approximately  equal.  They  consistently  applied  this  concept,  developing  expressions 
for  arrival  time  means  and  variances  for  nonreactive  solutes  in  the  "Standard  Model 
Aquifer,"  a  hypothetical  aquifer  characterized  by  a  heterogeneous  conductivity  field 
with  a  prescribed  correlation  structure  and  subject  to  certain  prescribed  boundary 
conditions  (see  Section  2.1).  These  equations  predicted  that  the  mean  arrival  time 
for  solute  subject  to  a  uniform  resident  injection  is  given  by 


where  angle  brackets  <>  denote  an  ensemble  average  operator  and  Vh  is  the  harmonic 
mean  Lagrangian  velocity.  From  the  small  deviation  assumption,  the  Lagrangian 
and  Eulerian  velocities  are  approximately  the  same,  thus  the  estimated  harmonic 
mean  Lagrangian  velocity  is  equal  to  the  estimated  harmonic  mean  Eulerian  velocity. 
However,  the  mean  arrival  time  for  uniform  resident  injection  of  solute  after  travelling 
several  integral  scales  is  given  by 


where  U  is  the  arithmetic  mean  Eulerian  velocity.  This  expression  is  also  the  "con- 
sistent first-order  approximation"  (CFOA)  of  the  travel  time.  Widespread  adoption 
of  the  CFOA  through  out  the  stochastic-Lagrangian  transport  literature  probably 
stems  Feynman's  requirement:  we  must  reproduce  what  we  already  know.  The  har- 
monic mean  of  a  positive  definite  process  is  always  less  than  the  arithmetic  mean,  so 
Equation  (1.1)  systematically  over-predicts  travel  times  for  large  displacements.  Dis- 
turbing was  that  the  estimators  derived  by  Shapiro  and  Cvetkovic  [1988]  using  clear, 
consistent,  intuitive  and  logical  methods  seemed  to  violate  Feynman's  rule,  especially 
in  light  of  the  success  enjoyed  by  the  Dagan  school. 


<r>=  x/Vh 


(1.1) 


<T>=  x/U 


(1.2) 
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This  harmonic/arithmetic  mean  discrepancy  was  explicitly  noted  by  Dagan  et  al. 
[1992],  who  relegated  the  Shapiro  and  Cvetkovic  [1988]  result  to  applicability  in  areas 
"close  enough  to  the  input  zone"  [Dagan  et  al,  1992  p.  1374]. 

Cvetkovic  et  al.  [1996]  developed  a  semi-analytical  expression  that  described  the 
transition  from  the  "near  field"  harmonic  mean  velocity  to  the  "far  field"  expression 
based  upon  the  these  known  "endpoints"  and  an  estimate  of  the  Lagrangian  velocity 
integral  scale.  Their  detailed  analysis  of  the  Lagrangian  velocity  field  revealed  that 
the  Lagrangian  velocity  is  nonstationary  in  displacements  in  space,  and  resulted  in  a 
nonlinear  propagation  of  the  mean  arrival  time  with  distance. 

This  was  strange.  Why  was  the  Lagrangian  velocity  field  nonstationary  when 
the  Eulerian  was  stationary?  It  was  commonly  assumed  that  the  stationarity  of  one 
implied  the  stationarity  of  the  other.  What  had  started  as  a  simple  preliminary 
literature  review  resulted  in  a  quandry.  The  paths  followed  in  pursuit  of  this  quandry 
led  to  the  "research  goals"  of  the  dissertation. 

1.3    Research  Goals 

There  were  four  broad  objectives  specified  for  this  work.  All  centered  around  the 
flow  of  water  and  the  transport  of  solutes  in  aquifers  characterized  by  heterogeneous 
hydraulic  conductivity  fields. 

1.3.1    Lagrangian  Velocity  Mean  and  Covariance 

A  primary  goal  of  this  research  was  to  develop  a  "Lagrangian  covariance  function" 
for  the  "Standard  Model  Aquifer".  This  was  extended  slightly  to  the  mean  and 
covariance  of  the  space-stationary  velocities  along  equal  flow-weight  streamlines  and 
time-stationary  velocities  along  equal  area-weight  streamlines.  Quantification  of  this 
covariance  function  greatly  simplifies  the  development  of  equations,  or  estimators, 
that  describe  the  movement  of  water  and  solutes  in  the  heterogeneous  velocity  fields 
associated  with  the  "Standard  Model  Aquifer." 
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1.3.2  Injection  Mode  Analysis 

The  stochastic-Lagrangian  framework  is  convenient  for  the  analysis  of  the  effects 
of  the  method  by  which  solute  is  introduced  into  the  flow  field,  or  injection  mode, 
upon  the  subeqent  transport.  These  effects  of  injection  mode  upon  nonreactive  solute 
transport  are  characterized  by  the  mean  and  variance  of  a  breakthrough  curve,  or 
mass  arrival  distribution,  and  by  the  mean  and  variance  of  mass  displacement.  A 
goal  of  this  research  was  to  develop  estimators  for  these  quantities. 

1.3.3  Heterogeneous  Reaction  Parameter 

Set  as  a  goal  were  estimators  for  the  mean  and  variance  of  a  heterogeneous  "reac- 
tion parameter"  integrated  along  trajectories  of  stationary  velocty.  These  estimators 
were  used  to  evaluate  the  coupled  effects  of  injection  mode  and  a  heterogeneous  linear 
equilbirium  sorption  process  in  a  stochastic-Lagrangian  framework. 

1.3.4  Application  to  "Reality" 

A  further  goal  was  to  develop  a  "traditional"  estimator  for  the  center  of  mass 
of  a  continuously  injected  solute  plume  subject  to  first-order  decay  and  to  evaluate 
data  from  a  natural  attenuation  experiment  and  numerical  experiments.  Estimator 
performance  was  evaluated  using  the  results  of  the  prevous  developments. 


CHAPTER  2 

LAGRANGIAN-EULERIAN  FLOW  FIELD  RELATIONSHIP 


Develop  an  intuitive  judgment  and  understanding  for  everything. 

Miyamoto  Musashi1 

2.1  Introduction 

Consider  the  steady,  well-developed  flow  of  water  in  an  aquifer  meeting  Bear's 
[1972]  definition  of  a  porous  medium.  In  most  natural  settings,  gentle  gradients 
induce  a  flow  well-quantified  by  Darcy's  law.  Knowledge  of  the  water  fluid  properties, 
relative  permeability  of  porous  medium,  and  potential  gradients  in  space  and  time 
would  assure  a  reasonable  estimate  corresponding  water  fluxes.2 

One  way  to  record  this  knowledge  would  be  to  create  maps  of  these  different 
properties  in  space  at  different  times.  This  framework  is  commonly  called  the  Eulerian 
framework,  and  is  the  most  common  approach  to  hydrologic  problems.  Another  way 
to  record  this  knowledge  is  to  map  properties  along  moving  coordinates  in  time. 
This  approach  is  commonly  called  the  Lagrangian  framework,  and  its  application 
may  without  controversy  be  described  as  less  widespread  than  Eulerian  approaches. 
Dagan  [1989]  and  Gelhar  [1993]  review  development  of  flow  of  fluids  and  transport  of 
contaminants  in  porous  media  in  both  of  these  frameworks. 


1  Musashi  [1982  p.  49] 

2  The  concepts  associated  with  the  following  discussion  are  quite  general,  and  are 
applicable  to  a  much  larger  variety  of  environmental  transport  problems.  The  em- 
phasis upon  groundwater  follows  from  an  interest  in  simplicity. 
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Consider  now  flow  of  water  in  a  more  specific  model  aquifer:  the  steady,  irrota- 
tional,  and  Darcian  flow  of  water  of  homogeneous  fluid  properties  in  a  finite  rectangu- 
lar two-dimensional  aquifer  of  heterogeneous  hydraulic  conductivity  (see  Figure  2.1). 
Mobile  water  9  completely  saturates  the  homogeneous  pore  space.  An  orthogonal 
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Figure  2.1:  Schematic  representation  of  an  aquifer. 
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Cartesian  coordinate  system  lies  parallel  to  the  aquifer  boundaries,  and  x  and  y  de- 
note its  directions.  Parallel  and  impervious  boundaries  lie  at  y  =  0  and  y  —  Ly. 
Constant  potential  boundaries  located  at  x  =  0  and  x  =  Lx  maintain  a  constant  av- 
erage potential  gradient  J  across  the  domain,  such  that  water  flows  into  the  aquifer  at 
x  —  0  and  out  of  the  aquifer  at  x  =  Lx  at  a  volumetric  discharge  rate  Q  (dimensions 
[L3T-1]).  For  an  intuitive  dimensional  consistency,  allow  the  aquifer  some  small  and 
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uniform  thickness  in  the  vertical  direction,  say  h.  The  total  area  of  the  vertical  plane 
at  x  =  0  through  which  water  enters  the  aquifer  is  A.  Define  the  spatially-averaged 
groundwater  seepage  velocity  at  this  inlet  plane  (IP)  as 

jj  I  A  ^a\da 

Lda  (2.1) 
=  Q_ 
OA 

For  a  large  class  of  hydrologic  applications,  water  quantity  is  the  sole  concern. 
For  these  problems,  the  "route"  taken  by  the  water  is  of  little  interest,  as  is,  say 
the  distribution  of  ages  or  residence  times.  Spatial  and  temporal  distributions  of 
hydraulic  properties  hold  some  interest,  but  only  in  how  they  relate  to  effective  bulk 
properties. 

However,  a  large  class  of  practical  hydrologic  problems  exist  for  which  the  trajec- 
tories of  water  parcels  across  the  flow  domain  is  of  great  importance.  Many  of  these 
problems  involve  the  fate  and  transport  of  contaminants  through  the  aquifer  system. 
The  fundamental  problem  is  this:  to  relate  an  Eulerian  flow  field  to  its  Lagrangian, 
or  trajectory-based,  equivalent. 

2.2    Flow  Field  Partitioning  Strategies 

A  streamline  is  a  hydrodynamic  entity  that  is  at  all  points  tangent  to  the  ve- 
locity vector.  For  two-dimensional  flow,  a  streamtube  is  an  object  defined  by  two 
streamlines.  Stream  surfaces,  the  intersection  of  which  form  streamlines,  define  three- 
dimensional  streamtubes.  Infinite  collections  of  streamlines  exist  that  subdivide  the 
flow  field  into  collections  of  streamtubes.  Again,  the  interest  of  clarity  suggests  a  focus 
upon  two-dimensional  flow.  For  a  lucid  discussion  of  streamlines,  stream  functions, 
and  stream  surfaces  in  two  and  three-dimensional  flows,  see  Bear  [1972]. 

Impervious  boundaries  are  streamlines,  thus  the  entire  model  domain  is  a  stream- 
tube  by  definition.  Viewed  as  a  "black  box,"  Eulerian  and  Lagrangian  distinctions 
have  little  meaning.  This  black  box  approach  is  the  purview  of  reactor  engineering, 
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an  important  and  well-developed  discipline,  from  which  a  great  deal  of  insight  and 
useful  mathematics  may  be  drawn.  For  an  excellent  review,  especially  with  respect 
to  concepts  related  to  transfer  functions,  see  Jury  and  Roth  [1990]. 

For  a  large  domain,  the  scale  at  which  chemical  reactions  takes  place  is  generally 
much  smaller  than  that  of  the  domain.  That  is,  processes  such  as  sorption,  microbial 
and  chemical  degradation,  and  the  like,  are  dependent  upon  chemical  potentials  across 
distances  on  the  order  of  the  pore  scale  to  perhaps  that  of  a  "micro-Darcy"  scale, 
defined  as  the  smallest  volume  for  which  Darcy's  law  is  generally  applicable.  However, 
the  scale  of  interest  tends  to  be  much  larger,  say  that  of  a  well,  or  its  capture  zone 
or  radius  of  influence,  or  perhaps  the  interface  across  which  groundwater  discharges 
into  a  surface  water  body,  or,  perhaps  the  Floridan  aquifer.  In  order  to  consider 
these  reactions  in  detail,  it  is  necessary  to  resolve  the  larger  domain  into  smaller 
sub-domains.  For  modeling  real  systems,  the  scale  and  degree  to  which  the  domain 
is  subdivided  should  be  based  upon  the  available  data. 

2.2.1    Simple  Subdivision 

Consider  the  water  flowing  into  the  model  aquifer  across  the  inlet,  or  injection, 
plane  (IP).  Assume  it  is  of  interest  to  divide  the  domain  into  two  sub-domains.  There 
are  many  ways  to  do  this.  One  would  be  to  divide  the  domain  "down  the  middle,"  that 
is,  at  a  line  parallel  to  the  x  axis  located  y  =  Ly/2  (curve  A  Figure  2.2).  This  division 
creates  two  sub-domains  of  equal  volume  and  gives  an  equal  "volume  weight"  to  each 
sub-domain.  The  division  also  creates  two  sub-IPs  of  equal  area  and  gives  an  equal 
"area  weight"  to  each  sub-domain.  For  a  perfectly  homogeneous  porous  medium,  the 
groundwater  velocity  vectors  will  be  parallel  to  this  line,  and  the  discharges  through 
each  sub-domain  will  be  the  same  as  well.  Moreover,  there  will  be  no  advective 
transfer  of  water  across  the  boundary.  However,  the  introduction  of  conductivity 
heterogeneity  in  the  transverse  direction  will  almost  certainly  lead  to  deviations  of 
the  streamlines,  and  to  streamlines  crossing  the  boundary.  This  implies  that  while 
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Figure  2.2:  Simple  subdivision  of  the  flow  field  by  one  of  four  strategies:  A)  equal- 
area  straight  line,  B)  equal-area  streamline,  C)  equal-flow  streamline,  D)  equal-flow 
straight  line. 

the  volume  of  water  contained  in  each  sub-domain  is  constant,  the  discharge  through 
a  plane  perpendicular  to  the  partition  is  a  function  of  plane  position  in  x. 

Consider  now  a  streamline  originating  at  (0,Ly/2)  (curve  B  on  Figure  2.2).  For 
a  homogeneous  medium,  this  corresponds  to  an  equal  volume  partitioning.  That  is, 
the  domain  centerline  is  a  streamline.  However,  transverse  heterogeneity  would  likely 
induce  changes  in  the  velocity  field  such  that  the  streamline  would  not  remain  on 
the  centerline.  For  general  transverse  heterogeneity,  the  flow  rate  entering  the  two 
domains  will  not  be  equal,  in  general.  Thus,  the  flow  weight  and  the  volume  weight 
of  the  two  sub-domains  are  unequal  in  this  case.  Again,  the  magnitude  of  the  areas 
through  which  water  enters  the  sub-domains  are  equal. 

Consider  now  a  streamline  originating  at  some  coordinate  (0,  y0)  such  that  the  vol- 
umetric flow  rate  is  equal  for  both  domains  (curve  C  on  Figure  2.2).  This  streamline 
corresponds  to  the  stream  function 

mv»]=mL>]:mo]  (2.2) 
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This  is  an  equal  partitioning  with  respect  to  the  flow  field,  or  the  stream  function. 
This  sort  of  partitioning  arises  quite  naturally  in  fluid  mechanics  and  in  numerical 
modeling  of  heterogeneous  flow  systems  using  streamtube  approaches  (see,  e.g.,  Theile 
et  al.  [1996]).  The  flow  weight  of  each  sub-domain  is  equal,  although  the  volume 
weights  are  different. 

A  fourth  scheme  would  be  to  extend  a  line  from  (0,y0)  to  (Lx,y0)  (curve  D  Fig- 
ure 2.2).  In  this  case,  the  sub-domains  share  only  equal  flow  rates  at  x  —  0.  In 
subsequent  discussions,  the  phrase  Eulerian  trajectory  will  refer  to  such  a  straight 
line  trajectory,  and  the  phrase  Lagrangian  trajectory  will  refer  to  trajectory  of  a 
streamline.  The  phrase  area  weight  will  refer  to  an  equal  area  partitioning  scheme 
at  a  definition  plane,  and  the  phrase  flow  weight  will  refer  to  an  equal  discharge  par- 
titioning scheme.  At  this  point,  it  is  important  to  note  that  no  mention  has  been 
made  of  boundary  conditions  as  they  relate  to  transport  equations.  These  partition- 
ing schemes  are  simply  ways  of  subdividing  a  flow  field.  An  area  weight  -  Eulerian 
trajectory  scheme  corresponds  to  traditional  spectral  approaches  that  assume  small 
deviations  about  a  mean  trajectory.  The  use  of  trajectory  is  quite  general,  and  does 
not  imply  a  "natural"  path  or  streamline. 

2.2.2    General  Subdivision 

Consider  dividing  the  flow  field  into  N  partitions  according  to  the  strategies  de- 
scribed (see  Figure  2.3).  Theile  et  al.  [1996]  demonstrated  that  streamlines  may 
represent  streamtubes.  Therefore,  the  remainder  of  this  discussion  will  focus  upon 
streamline  properties,  with  the  implication  that  these  properties  may  be  extended  to 
streamtubes  as  well.  A  primary  purpose  of  subdividing  the  domain  is  to  represent  the 
continuum  of  domain  properties  with  a  discrete  set  of  data.  A  regular  finite-difference 
grid  may  be  thought  of  as  a  set  of  points  at  equal  spacing  in  space  along  a  collection 
of  area-weighted  Eulerian  trajectories.  An  alternative  method  would  be  to  create 
spacings  based  upon  travel  times  along  the  trajectory,  where  travel  time  is  defined  as 
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B 

Figure  2.3:  Four  different  partitioning  strategies  for  the  same  flow  field:  A)  area 
weight  -  Eulerian  trajectory,  B)  area  weight  -  Lagrangian  trajectory,  C)  flow  weight 
-  Lagrangian  trajectory,  D)  flow  weight  -  Eulerian  trajectory. 
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the  path  integral  of  the  inverse  velocity  over  the  trajectory 


r[l]  =  / 
Jo 


ds 


(2.3) 


V[C]  •  a 


where  s  is  a  unit  vector  that  is  at  all  points  the  parallel  to  the  trajectory.  It  is  clear 
that  this  sort  of  Lagrangian  description  has  little  meaning  in  a  static  system,  and  is 
in  this  sense  less  general  than  Eulerian  descriptions.  Thus,  Lagrangian  description 
is  a  tool  reserved  for  dynamic  systems  and  is  a  supplement  to  Eulerian  descriptions. 
It  is  worthy  to  note  the  distinction  between  a  description,  a  general  approach,  and  a 
trajectory,  a  specific  entity  used  in  a  description. 

To  summarize,  a  primary  objective  is  to  approximate  a  steady  and  irrotational 
flow  field  continuum  with  a  discrete  set  of  trajectories.  The  trajectories  follow  either 
a  straight  line  (Eulerian  trajectory)  or  a  streamline  (Lagrangian  trajectory).  The 
trajectories  carry  either  an  equal  flow-weight  or  an  equal  injection  area-weight. 


Consider  once  again  the  model  aquifer  depicted  in  Figure  2.1.  As  mentioned  in  the 
introduction  of  this  chapter,  the  aquifer  has  certain  features  in  which  a  hydrologist 
might  be  interested,  e.g.,  the  intrinsic  permeability  or  the  porosity.  These  features 
may  be  thought  of  as  a  field  in  which  parameters  vary  continuously.  This  is  a  funda- 
mental concept  of  continuum  mechanics  and  is  discussed  in  a  hydrologic  context  by 
Bear  [1972],  Dagan  [1989],  and  Gelhar  [1993].  The  ultimate  description  of  any  such 
field  would  be  a  map  of  parameter  values  valid  at  any  point,  time,  or  scale.  In  lieu  of 
this  "ultimate  description,"  stochastic  hydrologists  characterize  spatial  heterogeneity 
with  conditional  probability  density  functions  that  reproduce  the  dominant  charac- 
teristics of  the  field  and  observed  values  at  given  locations.  The  entire  subject  of 
geostatistics  is  devoted  to  this  end  (see  Journel  and  Huijbregts  [1978]).  An  Eulerian 
field  refers  to  area-weighted  Eulerian-trajectory  oriented  field.  A  Lagrangian  field 
refers  to  a  flow- weighted  Lagrangian-trajectory  oriented  field. 


2.3    Descriptions  of  Fields 
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It  is  instructive  to  illustrate  these  concepts  with  an  example  and  to  interpret  the 
results  of  this  illustration. 


The  illustration  of  these  concepts  requires  a  heterogeneous  flow  field  and  col- 
lections of  trajectories,  both  Eulerian  and  Lagrangian,  that  correspond  to  this  flow 
field.  The  requirement  that  the  flow  field  be  based  upon  "reference  fields"  of  known 
properties  enhances  the  possibility  of  theoretical  interpretation  drawn  from  the  large 
body  of  literature  concerning  flow  and  solute  transport  in  idealized  porous  media.  A 
further  requirement  that  our  reference  fields  bear  some  resemblance  to  what  might 
be  expected  in  nature  enhances  the  possibility  of  useful  applicability. 

A  numerical  simulation  code  modeled  steady  groundwater  flow  in  a  hypothetical 
heterogeneous  two-dimensional  aquifer.  A  particle  tracking  code  traced  Eulerian  and 
Lagrangian  trajectories  across  this  aquifer,  recording  data  that  included  the  local 
groundwater  velocity  and  local  hydraulic  conductivity  at  intervals  in  space  and  time 
along  the  different  trajectories.  Statistics,  specifically  mean,  variance,  and  covariance, 
summarize  the  results  of  these  "measurements." 

The  conceptual  view  of  the  aquifer  is  that  given  by  Figure  2.1.  To  generalize  the 
system,  correlation  length,  or  integral  scale,  of  the  log  conductivity  field  Ainfc  serves 
as  a  characteristic  length  with  which  to  normalize  length  scales.  The  characteristic 
time  \\nk/U  serves  to  normalize  time  scales,  where  U  is  the  arithmetic  mean  velocity 
defined  in  Equation  (2.1).  The  definition  of  an  integral  scale  is  the  integral  from  zero 
to  infinity  of  the  correlation  function  R,  or, 


2.4    Numerical  Experiments 


(2.4) 
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The  design  criteria  for  the  model  aquifer  are  as  follows 
In  A:  correlation  length  X\nk  =  1 
effective  conductivity  Ke  =  1 
average  gradient  J  =  1  /50 
mobile  water  content  6  =  1/5 

From  these  criteria  and  Darcy's  law,  the  area-averaged  velocity  is  U  =  KeJ/6  =  1/10. 
Clearly,  the  parameter  values  are  artificial  and  selected  for  convenience,  but  values 
for  Ke  in  md~l  and  6  fall  in  ranges  typical  for  sand  [Domenico  and  Schwartz,  1990]. 

Two  5000  replicate  Monte  Carlo  simulations  of  aquifers  with  log  conductivity 
variances  of  1/2  and  1,  respectively,  comprised  the  suite  of  aquifer  simulations.  The 
following  sections  describe  the  experimental  design  and  the  results. 

2.4.1    The  Log  Conductivity  Field 

The  foundation  upon  which  to  build  the  illustration  is  a  heterogenous  hydraulic 
conductivity  field.  A  commonly-used  stationary  lognormally  distributed  and  expo- 
nentially correlated  isotropic  model  serves  as  the  Eulerian  hydraulic  conductivity 
reference  field  (see  e.g.,  Freeze  [1975]  or  Gelhar  and  Axness  [1983]).  Working  from  a 
hypothesized  conductivity  field  is  preferable  to  generating  a  hypothetical  velocity  field 
with  specified  characteristics  (e.g.,  Bellin  et  al.  [1994])  because  an  objective  of  this 
work  is  to  describe  the  underlying  conductivity  field  along  Lagrangian  trajectories 
and  relate  its  Lagrangian  description  to  its  Eulerian  properties. 

For  completeness,  it  is  pertinent  to  briefly  review  description  of  a  heterogeneous 
field  by  treating  the  process  as  being  "random"  Measured  values  of  hydraulic  con- 
ductivity of  virtually  any  natural  porous  formation  will  vary  from  place  to  place. 
In  the  absence  of  measurement  error,  this  variability  in  the  measured  values  reflects 
the  heterogeneity  of  myriad  physical  factors  and  processes  that  control  the  hydraulic 
conductivity.  As  mentioned,  measurements  of  hydraulic  conductivity  in  the  field  of- 
ten follow  a  lognormal  distribution.  In  these  cases,  it  is  useful  to  conceive  of  the 
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conductivity  field  as  a  realization  of  a  lognormal  random  process.  The  process  is 
said  to  be  stationary  if  the  ensemble  mean  and  variance  do  not  vary  in  space.  The 
process  is  nonstationary  if  some  trend  exists  in  the  mean  or  variance.  In  practice, 
only  single  realization  is  available,  and  the  "ensemble"  properties  are  ultimately  un- 
knowable. However,  if  the  field  is  large  compared  to  the  scale  of  variability,  one  may 
assume  ergodicity,  or,  that  the  ensemble  parameters  may  be  estimated  from  a  sample 
drawn  from  the  single  realization  since  there  exists  a  the  possibility  of  a  large  num- 
ber of  statistically-independent  observation  points.  For  a  more  detailed  development 
of  these  concepts  from  the  hydrologist's  perspective,  see  Gelhar  [1993].  For  a  more 
formal  development  of  these  concepts,  see  Papoulis  [1991]. 

The  model  of  hydraulic  conductivity  variability  is  a  stationary,  exponentially- 
correlated  isotropic  lognormal  distribution.  A  lognormal  random  variable  is  one  for 
which  its  natural  logarithm  is  normally  distributed.  As  such,  three  parameters  sum- 
marize the  statistical  properties  of  the  adopted  conductivity  model,  namely  the  mean 
of  log  k  /link,  the  variance  of  log  k  afnfc,  and  the  correlation  length  of  log  k  \\nk-  The 
covariance  model  is  given  by 


where  r  is  the  separation  between  two  points  in  the  log  conductivity  field.  Use  of 
this  simple  covariance  model  requires  stationarity  of  variance,  and  usually  implies  a 
stationarity  of  the  mean  as  well.  The  moments  of  the  "real  space"  conductivity  field 
are  given  by 


where  n  is  the  order  of  the  moment. 

In  the  absence  of  a  perfectly  described  physical  aquifer,  the  method  of  turning- 
bands  generated  equally-likely  realizations  of  a  hypothetical  conductivity  field  with 


C\nk[r]  =  af^expf-r/Ainfc] 


(2.5) 


<&n>=  exp[n//inJt  +  n2afnk/2] 


(2.6) 


16 

prescribed  ensemble  statistics  [Tompson  et  al,  1989].  Turning-bands  was  capable 
of  producing  two  and  three-dimensional  fields  with  the  desired  statistical  proper- 
ties and  was  available  as  a  portable  and  easily  modified  code.  Figure  2.4  illustrates 
the  lognormality  of  the  a  typical  two-dimensional  log  conductivity  field  realization. 
Figure  2.5  compares  an  estimate  of  the  longitudinal  covariance  to  the  hypothesized 
model.   In  order  to  capture  the  features  of  the  log  conductivity  spatial  variability,  the 


-3-2-10  1  2  3 

log  conductivity 

Figure  2.4:  Log  conductivity  distribution  from  a  41400  node  realization  generated 
by  the  turning-bands  method  compared  to  the  standard  normal  distribution.  Target 
distribution  parameters  are  fj,\nk  =  0  and  afnk  =  1. 

turning-bands  algorithm  generated  five  log  conductivity  observations  per  correlation 
length.  This  level  of  discretization  falls  in  the  range  of  other  similar  sets  of  simula- 
tions. Cvetkovic  et  al.  [1998]  used  a  discretization  of  4  nodes  per  Ainfc  for  simulations 
of  fields  with  ofnfc  as  large  as  1.  Bellin  et  al.  [1992]  used  a  discretization  of  8  nodes 
per  Ainfc  for  simulations  of  fields  with  a^nk  as  large  as  1.6. 
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Figure  2.5:  Estimated  log  conductivity  longitudinal  covariance  estimated  from  a  re- 
alization generated  by  the  turning-bands  method  compared  to  an  exponential  covari- 
ance function.  Variance  of  log  conductivity  afnk  —  1  and  correlation  length  \\^  =  1. 
Spacing  for  log  conductivity  values  was  Ainfc/5. 

Trajectories  covered  several  characteristic  lengths  X\nk  and  several  characteristic 
times  \\nk/U  to  capture  both  near-field  and  far- field  behavior  in  space  and  Lagrangian 
time.  The  first  requirement  was  easy  to  enforce  by  making  the  domain  size  equal  to  or 
greater  than  the  required  number  of  correlation  lengths.  The  second  requirement  was 
somewhat  more  problematic,  since  the  distance  traveled  in  some  number  of  \\nk/U 
varies  from  trajectory  to  trajectory.  Preliminary  simulations  indicated  that  a  domain 
40Ainfc  long  was  adequate  to  capture  travel  up  to  10  characteristic  time  scales  \\nk/U 
in  2000  trajectory  realizations  generated  in  the  most  heterogeneous  field  considered 
(ofnjt  —  1).  Thus,  the  selected  domain  extent  was  40Ainjt. 

The  conceptual  model  aquifer  is  a  large,  bounded,  primarily  two-dimensional, 
steady  flow  field.  As  such,  it  should  exhibit  a  variety  of  head,  head  gradient,  and 
velocity  nonstationarities  due  to  the  boundary  conditions  (see  Figure  2.1).  While 
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researchers  such  as  Rubin  and  Dagan  [1988], [1989],  and  Osnes  [1998]  have  investi- 
gated such  nonstationarities,  the  additional  complication  of  a  rigorous  consideration 
of  these  sorts  of  nonstationarities  would  unduly  obfuscate  the  very  simple  concepts 
illustrated  here.  To  minimize  the  effect  of  these  nonstationarities,  "solute  transport" 
is  constrained  to  take  place  some  distance  from  the  boundaries  of  the  numerical  do- 
main. Following  the  example  of  Bellin  et  al.  [1992],  a  3Ainfe  "buffer"  was  established 
on  both  sides  of  the  required  40  for  a  total  longitudinal  domain  length  of  46Ainfc. 
Preliminary  simulations  indicated  that  the  maximum  transverse  displacement  of  a 
streamline  in  a  afnk  =  1  field  was  around  one  third  the  longitudinal  displacement. 
For  a  streamline  originating  on  the  aquifer  centerline  parallel  to  the  x  axis,  an  aquifer 
width  of  36Ainfc  was  deemed  adequate  to  contain  large  transverse  displacements  over  a 
40Ainfc  longitudinal  displacement  and  provide  a  similar  "buffer"  to  mitigate  boundary 
effects. 

Combining  the  "physical"  domain  requirements  with  the  discretization  required 
to  capture  the  log  conductivity  variability  resulted  in  a  numerical  grid  230  by  180 
for  a  two-dimensional  aquifer  and  a  numerical  domain  with  41400  nodes.  A  three- 
dimensional  analysis  under  the  same  assumptions  would  have  required  a  numerical 
domain  of  more  than  7  million  nodes.  Simulations  of  this  magnitude  were  untenable 
in  light  of  the  available  computational  resources. 

2.4.2    Velocity  Field  Generation 

A  mixed  finite-element  scheme  solved  coupled  velocity  and  pressure  equations  de- 
rived from  Darcy's  law  and  the  continuity  equation.  The  code  took  input  parameters 
that  matched  the  design  criteria  and  the  log  conductivity  field  generated  by  turning- 
bands.  Andrew  I.  James  wrote  and  tested  the  core  of  the  code.  James  et  al.  [1997] 
used  a  derivative  version  of  this  code  to  generate  flow  fields  for  synthetic  tracer  tests. 
A  mixed  finite-element  solution  to  the  flow  field  was  specified,  following  the  results 
and  arguments  of  Mose  et  al.  [1994].   The  mixed  finite-element  scheme  is  known 
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to  provide  an  accurate  velocity  field  solution  for  groundwater  flow  in  heterogeneous 
conductivity  fields,  one  that  is  more  accurate  than  other,  more  traditional,  schemes 
[Mose  et  ai,  1994]. 

2.4.3  Trajectory  Generation 

An  objective  of  this  research  was  to  record  data,  such  as  local  velocity  and  and 
conductivity,  along  different  collections  of  trajectories,  both  Eulerian  and  Lagrangian. 
These  trajectories  were  traced  using  a  "particle-tracking"  scheme  based  upon  a  semi- 
analytical  method  presented  by  Pollock  [1988].  This  method  assumes  that  velocity 
varies  linearly  over  an  element,  and  that  the  velocity  in  some  direction,  say  i,  is 
independent  of  displacements  in  orthogonal  directions.  From  these  assumptions,  it  is 
a  relatively  simple  exercise  to  determine  the  displacement  of  a  "particle"  in  a  linear 
velocity  field  for  a  given  time,  or  the  time  required  to  make  a  certain  displacement. 
A  computer  code  recorded  particle  position,  time,  total  path  length  along  trajectory, 
velocity,  local  conductivity,  and  a  value  from  another  random  field  with  the  same 
statistical  properties  as,  but  uncorrelated  to,  the  log  conductivity  field  at  several 
control  planes  equally-spaced  along  the  mean  flow  direction  and  at  several  "control 
times"  after  injection  along  both  trajectories.  This  uncorrelated  random  variable  was 
recorded  for  reactive  transport  simulations,  to  be  discussed  in  subseqent  chapters. 

2.4.4  Monte  Carlo  Simulation 

The  possibility  of  theoretical  interpretation  of  the  results  is  enhanced  by  requiring 
near-ergodic  samples  from  the  population.  That  is,  the  samples  must  reproduce  the 
population  statistics.  To  do  this,  the  number  of  statistically  independent  replicates 
drawn  from  the  velocity  and  conductivity  fields  must  be  quite  large.  Rather  than 
attempt  a  massive  single-realization,  the  following  method  generates  a  large  number 
of  statistically  independent  realizations 
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•  generate  a  large,  two-dimensional  conductivity  field  with  the  desired  statistical 
properties 

•  generate  a  flow  field  by  imposing  a  constant  mean  gradient  across  the  saturated 
conductivity  field 

•  trace  the  trajectories  of  two  particles  injected  at  a  common  "control  location", 
one  trajectory  is  a  straight  line  or  Eulerian  trajectory  and  the  other  a  streamline 
or  Lagrangian  trajectory,  and  record  properties  at  equal  lags  in  space  and  in 
time  along  trajectories  that  pass  through  "inner  domains"  that  are  relatively 
undisturbed  by  nonstationarities  induced  by  the  boundaries 

At  this  juncture,  it  is  perhaps  worth  mentioning  a  minor  "philosophical"  point. 
These  simulations  do  not  represent  a  real  aquifer,  per  se,  but  rather  a  hypothetical 
aquifer  with  known  properties.  The  goal  is  understand  the  convolution  of  the  physics 
of  the  groundwater  flow  and  a  known  heterogeneity  that  has  analogs  in  nature,  and  it 
is  hoped  that  the  reader  finds  some  merit  to  this  modest  goal.  The  41400  conductivity 
values  needed  for  any  single  realization  do  not  represent  41400  measurements  taken 
from  some  real  aquifer.  The  Monte  Carlo  aspect  of  this  experiment  reenforces  this 
point. 

2.4.5    Sampling  Strategy 

Consider  "dropping  particles"  into  a  flow  field,  with  no  particular  regard  to  lo- 
cation. These  particles  move  along  streamlines  without  interaction  with  the  porous 
medium.  Any  one  location  into  which  one  is  dropped  is  as  equally-likely  as  any  other. 
The  trajectories  associated  with  these  particles  carry  an  equal  area  weight,  analogous 
to  equal-area  subdivision  of  an  injection  plane.  Consider  now,  the  same  disposal  of 
particles,  but  this  time  an  area  is  assigned  to  each  particle  at  the  injection  point  such 
that  an  equal  discharge  is  associated  with  each  particle,  and  its  associated  trajectory. 
For  a  homogeneous  mobile  water  content,  this  area  is  inversely  proportional  to  the  lo- 
cal velocity  normal  to  the  area.  Thus,  each  particle  may  be  thought  to  carry  an  equal 
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flow  weight.  In  this  way,  it  is  possible  to  "map"  from  one  collection  of  trajectories  to 
another. 

At  a  given  control  plane  or  time,  the  area-weighted  mean  and  variance  of  some 
property  over  N  realizations,  say  v  was  calculated  using  the  familiar  estimators3 

1  N 

i=i 

1  N 

vark]  =  J^(Y1  vi  ~  N^  (2-8) 
i=i 

While  these  particular  estimators  seem  obvious,  it  is  important  to  realize  that  "drop- 
ping the  particle"  into  the  flow  field  assigns  an  equal  area  weight  to  each  trajectory. 
For  a  single  realization,  this  is  equivalent  to  starting  each  particle  with  equidis- 
tant spacing.  These  estimators  are  based  upon  an  assumption  of  equally-likely 
statistically-independent  samples  (see,  e.g.,  Mood  et  al.  [1974]).  In  order  to  draw 
even  a  few  hundred  independent  samples  from  a  single  realization,  the  required  flow 
domain  would  be  enormous.  This  was  a  primary  consideration  for  adopting  Monte 
Carlo  simulation  using  a  minimal  domain. 

The  flow-weighted  mean  property  v  was  calculated  via 

Ei=i«<[o] 


3  Press  et  al.  [1988]  have  criticized  the  variance  estimator  in  Equation  (2.8)  as 
being  prone  to  accumulation  of  round-off  error  for  large  data  sets.  The  rather  sub- 
stantial size  of  the  data  sets  are  naturally  a  concern.  Comparision  of  the  variance 
calculated  by  a  corrected  two-pass  method  presented  by  Press  et  al.  [1988]  and  the 
"traditional"  estimator  yielded  no  difference  between  the  two  methods  within  the 
first  five  significant  digits. 
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where  u;[0]  is  the  initial  velocity  recorded  for  trajectory  i.  The  variance  by 

Ei=i  w 

An  alternative  method  for  generating  a  collection  of  streamlines  would  be  to  iden- 
tify the  streamline  corresponding  to  the  "center"  stream  function  expressed  in  Equa- 
tion (2.2).  In  this  case,  the  estimators  of  the  equal  flow  statistics  would  be  Equa- 
tions (2.7)  and  (2.8).  The  equal  area  statistics  would  be  given  by 

_  =E£Wm[o]  (211) 
Et=i  i/«i[o] 

and 

varM  =  ^tM-^lum 
E£i.l/«*[0] 

The  "reference"  collection  of  the  trajectories  selected  uniformly  in  space  was  specified 
in  the  simulations  for  convenience. 

2.4.6    Simulation  Results 

The  following  sections  present  the  results  of  two  5000  replicate  Monte  Carlo  simu- 
lations. Detailed  discussion  of  these  results  is  deferred  until  after  the  development  of 
some  theoretical  models  with  which  to  analyze  these  results.  However,  some  heuris- 
tics are  provided  to  familiarize  the  reader  with  some  of  the  concepts  developed  in 
subsequent  sections  and  chapters. 

A  summary  of  the  notation  used  with  the  following  results  provides  clarity.  AW 
denotes  an  area- weighted  trajectory  and  FW  denotes  a  flow- weighted  trajectory  .  ET 
denotes  an  Eulerian  trajectory  and  LT  denotes  a  Lagrangian.  A  particular  weight- 
trajectory  combination  is  concatenated  the  with  a  /.  Thus,  an  area-weighted  Eulerian 
trajectory  is  denoted  as  AW/LT. 

For  a  given  simulation,  that  is,  for  one  set  of  5000  replicates  for  a  log  conductivity 
field  of  some  given  a?nk,  the  statistics  of  the  log  conductivity  and  velocity  along  the 
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four  different  trajectories  are  calculated  at  intervals  in  time  and  in  space.  Thus,  for 
any  one  property,  say  log  conductivity,  and  any  one  statistic,  say  the  mean,  there  will 
be  two  sets  of  results,  namely  that  for  afnk  =  1/2  and  that  for  afnk  =  1.  Associated 
with  each  set  will  be  two  subsets,  namely  statistics  along  the  four  trajectories  recorded 
at  space  or  time  intervals. 

For  convenience,  the  notation  and  the  trajectories  to  which  each  notation  corre- 
sponds as  depicted  in  Figures  2.2  and  2.3  is  summarized  here: 
AW/ET  area- weighted  Eulerian  trajectory  (Figures  2.2  and  2.3  A) 
AW/LT  area- weighted  Lagrangian  trajectory  (Figures  2.2  and  2.3  B) 
FW/LT  flow-weighted  Lagrangian  trajectory  (Figures  2.2  and  2.3  C) 
FW/ET  flow- weighted  Eulerian  trajectory  (Figures  2.2  and  2.3  D) 

2.4.7    Log  Conductivity  Statistics 

The  log  conductivity  field  serves  as  the  "foundation"  upon  which  the  simulations 
are  built,  as  well  as  the  theory.  The  data  for  the  area- weighted  Eulerian  trajectories 
correspond  to  those  design  criteria  specified.  The  log  conductivity  mean  and  vari- 
ance for  AW/ET  appear  stationary  in  space  (see  Figures  2.6  and  2.7).  That  is,  they 
appear  to  "bounce  around"  their  "target"  values  with  no  apparent  trend.  Moreover, 
the  statistics  of  FW/LT  appear  likewise  stationary,  but  the  mean  is  somewhat  larger 
and  the  variance  slightly  smaller  than  AW/ET.  The  magnitude  of  the  differences 
appear  proportional  to  ofnfc.  A  comparison  of  the  "target"  values,  e.g.,  the  specified 
inputs  into  the  turning-bands  algorithm,  and  the  observed  ensemble  values  along  the 
trajectories  is  presented  in  Table  2.1.  The  "observed"  entry  in  Table  2.1  represents 
the  statistic  average  along  the  trajectory,  as  estimated  by  regressing  a  zero-slope  line 
through  the  data  using  a  generalized  least-squares  method.  There  is  close  correspon- 
dence between  the  target  and  observed  values,  supporting  an  ergodic  hypothesis  for 
the  log  conductivity  statistics.  That  is,  sample  statistics  are  an  accurately  estimation 


24 


0.6 


c 
=1 


0.5 
0.4 
0.3 
0.2 
0.1 
0 

-0.1 
-0.2 

0.6 
0.5 
0.4 
0.3 
0.2 
0.1 
0 

-0.1 
-0.2 


+__+  .  *  * 


*  » 


FW/LT  + 
AW/LT  x 
FW/ET  « 
AW/ET  □ 

Eq.  2.72   

reference   


+*  X       *       ix,  S  + 


+  +  +    **    ¥      5*  ** 

+     *  *  * 

x   x  * 
*    X  x 


I  L. .I*...-. JW.!  A. ..S.  "»fl-|J.l  

□  X 


FW/LT 
AW/LT 
FW/ET 
AW/ET 
Eq.  2.72 
reference 


10 


20 


30 


40 


Figure  2.6:  Mean  log  conductivity  as  function  of  displacement  in  mean  flow  direction 
for  the  four  different  trajectory  collections  and  afnk  =  1/2  (a)  and  afak  =  1  (b). 
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Figure  2.7:  Log  conductivity  variance  as  function  of  displacement  in  mean  flow  direc- 
tion for  the  four  different  trajectory  collections  and  afnk  =  1/2  (a)  and  a?ak  =  1  (b). 
Variances  are  normalized  by  target  variance  of  Eulerian  log  conductivity  field  ofnk. 
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Table  2.1:  Comparison  of  target  log  conductivity  values  and  observed  values  along 
area-weighted  Eulerian  trajectories. 


of  the  population  statistics.  Again,  this  is  important  for  a  quantitative  theoretical 
analysis  of  these  results. 

The  AW/LT  and  the  FW/ET  exhibit  distinctly  nonstationary  behavior  in  the 
mean  and  less  distinct  behavior  in  variance  as  a  function  of  displacement  along  the 
mean  flow  direction  (see  Figures  2.6  and  2.7).  At  the  injection  point,  the  AW/ET  and 
AW/LT  statistics  are  identical  and  the  FW/ET  and  FW/LT  statistics  are  identical, 
because  these  trajectories  carry  the  same  weight  and  the  points  in  consideration  are 
identical  (see  e.g.,  Figure  2.2).  At  very  large  displacements,  one  would  expect  the 
statistics  of  like  trajectories  to  assume  similar  characteristics,  regardless  of  weight. 
As  the  trajectory  traverses  several  integral  scales  that  characterize  the  heterogeneity 
in  question,  "information"  about  its  starting  location,  is  diminished.  In  the  case  of 
the  FW/ET,  the  statistics  "start"  at  the  same  place  as  FW/LT,  and  eventually  "end" 
at  the  same  place  as  the  AW/ET.  Similarly,  the  AW/LT  statistics  start  where  do  the 
AW/ET  and  end  where  do  the  FW/LT.  In  fact,  the  data  exhibit  this  behavior  (see 
Figures  2.6  and  2.7). 

The  log  conductivity  statistics  in  time  show  somewhat  different  behavior.  Only 
the  AW/LT  appear  stationary  in  time  (see  Figures  2.8  and  2.9).  There  is  an  impor- 
tant observation  to  be  gleaned  from  these  results.  The  AW/LT  statistics,  in  addition 
to  appearing  stationary,  also  appear  to  have  a  similar  value  to  the  input  AW/ET 
statistics  (see  Table  2.2).  Recall,  the  AW/LT  statistics  are  clearly  nonstationary  in 
space,  and  are  clearly  distinct  from  those  of  the  AW/ET.  This  observation  greatly 
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Figure  2.8:  Mean  log  conductivity  as  function  of  travel  time  to  control  planes  along 
mean  flow  direction  for  the  four  different  trajectory  collections  and  afnk  =  1/2  (a) 
and  alk  =  1  (b). 
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Figure  2.9:  Log  conductivity  variance  as  function  of  travel  time  to  control  planes 
along  mean  flow  direction  for  the  four  different  trajectory  collections  and  afnk  =  1/2 
(a)  and  afnk  =  1  (b).  Variances  are  normalized  by  target  variance  of  Eulerian  log 
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Table  2.2:  Comparison  of  target  log  conductivity  values  for  area- weighted  Eulerian 
trajectories  in  space  and  observed  area-weighted  Lagrangian  trajectories  in  time. 

parameter  target  observed 


MinfcK^  =  1/2)  0  -0.000 
HinkKk  =  1)  0  -0.017 

1/2  0.497 
1  1.000 

simplifies  relating  the  Eulerian  and  Lagrangian  fields,  as  shall  be  shown.  The  "long- 
time" asymptotic  mean  log  conductivity  values  for  both  Eulerian  and  Lagrangian 
trajectories  are  somewhat  less  than  their  "large-displacement"  counterparts,  and  the 
magnitude  of  the  difference  is  proportional  to  o^nk  (compare  Figures  2.6  and  2.8). 
This  is  consistent  with  the  harmonic  averaging  implicit  with  a  time-averaged  time- 
dependent  velocity  and  the  arithmetic  averaging  implicit  with  a  space-averaged  space- 
dependent  velocity.  Recall,  we  map  the  quantities  onto  the  trajectories  by  following 
velocities,  whether  along  "natural"  streamlines  (Lagrangian  trajectories)  or  not  (Eu- 
lerian trajectories).  A  simple  heuristic  is  as  follows.  Low  conductivities  correspond  to 
low  velocities.  For  an  equal-time  spacing  (or  equally-likely  observations  in  time)  more 
"measurements"  will  come  from  low  velocities,  since  the  particle  tracing  the  trajec- 
tory spends  more  time  going  slow.  For  an  equal-distance  spacing  more  measurements 
will  come  from  high  velocities,  since  the  particle  "covers  more  ground"  while  going 
fast.  This  is  perhaps  more  intuitive  when  discussed  in  terms  of  the  velocity  field  itself. 

2.4.8    Velocity  Statistics 

The  mean  velocities  along  the  different  trajectories  exhibit  a  similar  behavior  as  do 
the  mean  log  conductivities  along  the  same  trajectories  (see  Figure  2.10  and  compare 
Figure  2.6).  The  AW/ET  exhibits  apparently  stationary  behavior  in  the  mean  and 
variance  (see  Figure  2.11)  in  accordance  to  established  theory  (e.g.,  Dagan  [1989]). 
The  AW/LT  exhibits  nonstationary  behavior  in  the  mean  as  observed  by  Cvetkovic 
et  al.  [1996].    Again,  the  AW/ET  and  FW/LT  trajectories  appear  to  be  stationary 
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Figure  2.10:  Comparison  of  position-dependent  mean  velocity  estimators  to  data  for 
the  four  collections  of  trajectories  for  afnk  =  1/2  (a)  and  afnk  —  1  (b).  Velocities  are 
normalized  by  U. 
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Figure  2.11:  Comparison  of  the  position-dependent  velocity  variance  estimators  to 
data  for  the  four  collections  of  trajectories  for  cr,2nfc  =  1/2  (a)  and  afnk  —  1  (b).  The 
variances  are  normalized  by  U2bafnk. 
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in  both  the  mean  and  variance.  The  velocity  variances  exhibit  behavior  to  similar 
to  that  of  the  means  for  their  corresponding  trajectories  and  quite  different  than  the 
behavior  exhibited  by  the  log  conductivity  variances  (compare  Figures  2.11  and  2.10 
and  compare  Figures  2.11  and  2.7).  There  is  a  qualitative  similarity  between  the 
nonstationary  velocity  means  and  the  nonstationary  velocity  variances.  A  heuristic 
is  that  the  velocity  statistic  "starts"  at  one  value  "in  the  near-field"  and  transitions 
to  an  asymptotic  value  "in  the  far-field."  The  following  sections  demonstrate  that 
these  endpoints  and  the  transition  are  predictable. 

The  qualitative  dissimilarity  between  the  log  conductivity  variance  and  velocity 
variance  is  an  artifact  of  examining  the  log  conductivity  values  on  the  one  hand  and 
the  "real  space"  velocity  values  on  the  other.  Consider  the  F W/LT  log  conductivity 
mean  and  variance  for  the  afnk  =  1  simulation.  The  FW/LT  mean  log  conductivity 
is  larger  than  that  of  the  AW/ET,  but  the  variances  are  roughly  the  same  value. 
In  fact,  regressing  a  zero-slope  line  through  the  data  results  in  a  "mean"  FW/LT 
variance  that  is  slightly  less  than  that  of  the  AW/ET.  However,  the  "real  space" 
conductivity  mean  and  variance  are  both  functions  of  both  the  mean  and  variance  of 
the  log  conductivity  process  (recall  Equation  (2.6)).  A  cursory  examination  of  the 
data  in  conjunction  with  Equation  (2.6)  reveal  that  the  mean  and  variance  of  the  "real 
space"  conductivity  along  FW/LT  are  both  higher  than  their  AW/ET  counterparts, 
and  the  magnitude  of  this  difference  is  proportional  to  afnk. 

The  mean  velocities  in  Figure  2.10  are  normalized  by  the  arithmetic  mean  veloc- 
ity U  predicted  from  the  input  parameters  for  the  flow  simulator.  That  the  mean 
velocity  along  AW/ET  closely  follows  1  is  a  good  indication  that  the  effective  con- 
ductivity expression  is  accurate.  Notice  that  the  FW/LT  and  far-field  AW/LT  values 
are  higher  than  AW/ET,  and  the  magnitude  is  proportional  to  the  variance  of  the 
log  conductivity  (see  Figure  2.10).  This  behavior  is  intuitively  correct.  Consider  a 
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particle  injected  into  a  flow  field.  Constraint  of  the  trajectory  to  an  Eulerian  tra- 
jectory forces  the  trajectory  through  low-velocity  areas  that  a  natural  streamline,  or 
Lagrangian  trajectory,  would  have  otherwise  bypassed.  Of  course,  some  streamlines 
pass  through  even  the  lowest  velocity  areas,  but  these  are  few  in  an  equal-flow  sense. 
Thus,  equal-area  weighting  preferentially  weights  low  velocity  zones.  That  there  are 
implications  for  sampling  schemes  in  heterogeneous  media,  especially  in  context  of 
multi-level  samplers,  should  be  intuitively  obvious  to  even  the  most  casual  observer. 
A  discussion  of  some  of  these  implications  is  made  in  a  subsequent  section,  after  these 
concepts  are  developed  further  in  a  more  quantitative  and  theoretical  framework. 

It  is  appropriate  to  evaluate  the  adopted  expression  for  the  AW/ET  longitudinal 
covariance  function  suggested  by  Dagan  and  Cvetkovic  [1993],  namely 

Cua[r)  =  alkb[e]U2  exp  [-b[e)r  / \Xnk]  (2.13) 

where  r  is  the  separation  between  points  along  the  ET,  and  b[e]  is  a  function  of  the 
anisotropy  ratio  e.  The  anisotropy  ratio  itself  is  a  function  of  the  correlation  lengths 
of  the  log  conductivity  in  the  different  directions.  The  primarily  consideration  is 
that  b  may  be  derived  from  theoretical  considerations  and  is  not  some  ad  hoc  fitting 
parameter.  The  value  of  this  function  for  a  two-dimensional  isotropic  aquifer  is  3/8. 
Again,  for  further  background  on  this  function,  see  Dagan  and  Cvetkovic  [1993].  In 
the  limit  r  =  0,  the  covariance  reduces  to  the  point  variance,  which  for  the  model 
case  is  a\  =  3[/2ofnfc/8.  For  the  remainder  of  this  discussion,  a  reference  to  b  contains 
its  dependence  upon  e  implicitly  The  subscript  ua  may  be  decoded  as  follows:  for 
a  subscript  tw,  the  t  is  the  trajectory  and  property  (u  for  Eulerian  velocity,  and  v 
for  Lagrangian  velocity)  and  the  w  is  the  weight  (a  for  area  and  /  for  flow).  Thus, 
Cua  is  the  covariance  function  for  velocities  along  Eulerian  trajectories.  We  present 
this  covariance  function  and  that  covariance  function  estimated  from  the  assumed- 
stationary  collection  of  AW/ETs  in  Figure  2.12  The  perceptive  observer  will  note 
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Figure  2.12:  Comparison  of  spatial  Eulerian  (AW/ET)  and  Lagrangian  trajectory 
(FW/LT)  velocity  covariance  data  to  covariance  model.  Data  are  taken  from  the 
ofnfc  =  1  set  of  realizations.  Covariances  presented  as  dimensionless  lag  correlation 
functions. 

that  the  analytical  expression  over-estimates  the  covariance  at  short  lags  and  under- 
estimates the  covariance  at  longer  lags.  This  discrepancy  and  its  implications  are 
discussed  in  some  detail  in  a  subsequent  section. 

The  velocity  variances  in  Figure  2.11  are  normalized  by  the  zero-lag  covariance  pre- 
dicted by  placing  the  input  simulation  parameters  into  Equation  2.13.  The  AW/ET 
velocity  variance  appears  well-predicted  by  this  expression,  as  evidenced  by  a  rela- 
tively good  fit  to  1  for  both  values  of  o?nk  in  Figure  2.11.  A  regression  of  a  zero-slope 
line  through  these  data  indicates  that  the  a  priori  estimate  under  predicts  the  ob- 
served values  by  about  five  percent  for  the  larger  o^nk  (see  Table  2.3)  The  apparent 
stationarity  of  the  velocities  along  the  FW/LTs  suggest  estimation  of  the  covariance 
function  for  these  velocities  as  for  the  AW/ET  (see  Figure  2.12).  The  correlation 
length  of  the  FW/LT  is  very  close  to  that  of  the  AW/ET,  although  the  actual  value 
is  probably  slightly  higher.  Recall  the  definition  of  the  correlation  length  as  given  by 
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Table  2.3:  Comparison  of  a  priori  Eulerian  velocity  variance  to  observed  values.  The 
decimal  representation  of  the  exact  fractional  values  of  the  predicted  variance  are 
provided  for  convenience. 

parameter  predicted  observed 

o2u{o'ik  =  1/2)  3/1600=0.001875  0.001871 
^(^i2nfe  =  1)  3/800=0.00375  0.00394 

Equation  (2.4).  The  correlation  function  at  each  considered  lag  along  the  estimated 
functions  is  systematically  larger  for  F W/LT,  even  though  the  magnitude  of  the  dif- 
ference is  small.  Integrating  these  functions  would  result  in  a  slightly  larger  value 
for  FW/LT,  and  thus  a  larger  correlation  length.  In  the  case  there  is  no  variability 
in  the  log  conductivity,  Eulerian  and  Lagrangian  trajectories  coincide.  Moreover,  an 
equal-area  weight  corresponds  to  an  equal  flow  weight.  Thus,  the  estimated  correla- 
tion functions  from  "either  set"  would  be  identical.  From  this  conclude  that  there 
is  some  dependence  of  the  correlation  length  of  the  FW/LT  trajectories  upon  afnk. 
Cvetkovic  et  al.  [1996]  observed  similar  behavior  in  log  velocity  correlation  functions 
estimated  from  AW/LTs  with  a  stationarity  assumption  for  the  nonstationary  AW/LT 
velocities.  Thus,  this  is  an  slight  extention  of  this  work  by  identifying  the  station- 
ary streamline-based  velocity  field  from  which  to  work,  and  work  from  that  field. 
Moreover,  this  supports  their  findings  of  an  increasing  correlation  length  with  afnk. 
Understanding  the  relationship  between  the  variability  of  the  log  conductivity  field 
and  the  correlation  length  of  the  associated  velocity  field  is  of  fundamental  impor- 
tance to  completely  understanding  the  nature  of  flow  in  heterogeneous  media.  Sadly, 
these  "secrets"  remain  uncovered  and  must  be  left  this  to  future  research  efforts. 

The  velocity  means  as  a  function  of  time  exhibit  qualitatively  similar  behavior 
as  do  the  log  conductivity  means  in  time  (see  Figure  2.13  and  compare  Figure  2.8). 
Again,  only  the  AW/LT  exhibits  stationary  behavior  in  time  for  both  the  mean  and 
variance  (see  Figure  2.14)  as  did  the  AW/LT  log  conductivity  values.  The  mean 
and  variance  of  the  AW/LT  in  time  are  approximately  those  of  the  AW/ET  in  space 
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Figure  2.13:  Comparison  of  time-dependent  mean  velocity  estimators  to  data  for  the 
four  collections  of  trajectories  for  afnk  =  1/2  (a)  and  a^nk  =  1  (b).  Velocities  are 
normalized  by  U. 
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Figure  2.14:  Comparison  of  time-dependent  velocity  variance  estimators  to  data  for 
the  four  collections  of  trajectories  for  afnk  —  1/2  (a)  and  afnk  =  1  (b).  The  variances 
are  normalized  by  U2bafnk. 


38 

for  both  simulations.  This  result  is  important,  because  it  implies  the  possibility  of 
mapping  from  AW/ET  to  AW/LT.  We  estimated  the  covariance  function  for  the 
stationary  AW/LT  velocities  (see  Figure  2.15).   As  an  estimate  of  the  covariance 
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Figure  2.15:  Comparison  of  temporal  Lagrangian  trajectory  (AW/LT)  velocity  covari- 
ance data  to  covariance  model.  Data  are  taken  from  the  afnfc  =  1  set  of  realizations. 
Covariances  presented  as  dimensionless  lag  correlation  functions. 


function,  the  substitution  r  =  Ut  was  made  in  Equation  2.13  and  the  result  compared 
to  the  estimated  covariance  function  (see  Figure  2.15).  The  closed-form  predictor 
behaves  similarly  to  its  spatial  counterpart:  early-time  correlation  is  under-predicted 
and  long-time  correlation  is  over-predicted  (compare  Figures  2.15  and  2.12). 

For  completeness,  one  other  aspect  of  the  velocity  field  was  examined,  namely  the 
distribution  of  initial  velocities  (see  Figure  2.16).  The  logarithm  of  the  initial  velocity 
was  rank-ordered  and  subject  to  the  transformation 

In  v  —  In  v 


\nv' 


J\n  v 


(2.14) 
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Figure  2.16:  Normalized  log  initial  velocities  for  two-dimensional  model  aquifer  com- 
pared to  the  standard  normal  distribution.  The  random  variable  was  normalized  by 
the  sample  mean  and  standard  deviation  to  illustrate  the  lognormality. 

where  lni;  is  the  sample  mean  and  Ein„  is  the  sample  standard  deviation.  Work- 
ing here  with  the  logarithm,  rather  than  the  velocity  itself,  better  demonstrates  the 
apparent  lognormality  of  the  velocity  fields.  Recall,  the  logarithm  of  a  lognormal 
random  variable  is  itself  normally  distributed.  Figure  2.16  indicates  that  the  trans- 
formed log  initial  velocities  are  well-described  by  a  standard  normal  distribution.  The 
lognormality  of  the  velocities  seems  logical  since  the  conductivity  field  is  itself  lognor- 
mally  distributed.  However,  the  velocity  is  in  general  a  product  of  random  variables: 
the  conductivity,  the  head  gradient,  and  the  mobile  water  content.  The  logarithm  of 
the  velocity  is  the  sum  of  random  variables.  Sums  of  random  variables  tend  toward 
normality  by  the  central  limit  theorem  (see  Mood  et  al.  [1974]),  and  thus,  products 
of  random  variables  tend  towards  lognormality.  It  would  be  interesting  to  test  this 
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hypothesis  by  examining  the  velocity  characteristics  of  flow  fields  arising  from  differ- 
ent log  conductivity  distributions  (e.g.,  fractal,  uniform,  white  noise,  etc.),  and  this 
is  left  to  future  work. 

2.4.9    Head  and  Head  Gradient 

Head  and  head  gradient  information  along  the  trajectories  was  not  recorded.  This 
is  a  regrettable  oversight,  since  much  of  the  behavior  of  the  flow  water  and  the  trans- 
port of  solutes  in  heterogeneous  media  can  be  summarized  as  the  interplay  of  head, 
head  gradient  and  conductivity  covariance  functions.  The  literature  concerning  Eu- 
lerian  head  and  head  gradient  fields  is  extensive  and  is  summarized  in  large  part  by 
Gelhar  [1993].  In  fact,  much  is  drawn  from  this  body  of  literature  for  some  of  the 
theoretical  developments  in  the  following  sections.  That  this  information  was  not 
recorded  detracts  mainly  from  the  completeness  of  this  work,  and  further  investiga- 
tions are  left  to  future  works. 

2.5  Theory 

2.5.1    Effective  Conductivity 

As  with  the  simulations,  the  theoretical  developments  are  founded  on  the  bedrock 
of  the  log  conductivity  field.  Section  2.4.7  demonstrates  that  the  simulation  results 
well  match  the  "target"  ensemble  quantities.  From  this  field  is  had  directly  a  derived 
parameter,  namely  the  a  priori  effective  conductivity  Ke.  An  expression  for  the 
effective  conductivity,  similar  to  that  given  by  Gelhar  and  Axness  [1983],  is  derived 
here  to  illustrate  the  mechanics  of  different  techniques  employed  through  out  these 
derivations.  Consider  a  one-dimensional  Darcy's  law  of  the  form 


(2.15) 


where  q  is  the  specific  discharge,  k  is  the  hydraulic  conductivity,  and  <j>  is  the  hydraulic 
head.  Consider  a  large  stationary  conductivity  field  such  that  the  log  conductivity 
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may  be  resolved  into  a  mean,  and  a  zero-mean  perturbation,  e.g.,  In  A;  =  fj,\nk  +  /• 
Let  the  head  gradient  likewise  be  resolved  into  a  mean  and  zero-mean  perturbation, 
e.g.,  -d(p/dx  =  J(l  +  ;').  Inserting  these  expressions  into  Equation  (2.15)  yields 

q  =  exp[/i[n  k  +  f]J{l  +  j) 

(2.16) 

=  KgJexp[f](l+j) 

Expand  the  exponential  term  in  this  equation  as  a  Taylor  series  about  zero,  truncating 
the  series  at  second-order,  and  expand  the  product  of  the  perturbation  terms 

q~KgJ(l  +  f  +  f/2)(l+j) 

(2.17) 

=  KgJ(l  +  f  +  f2/2  +  j+fj  +  f2j/2) 
Take  the  expected  value  of  this  expression,  dropping  terms  higher  than  second-order, 
resulting  in  an  expression 

<  q  >=  KgJ(l  +  al  k/2+  <  fj  >)  (2.18) 

Consider  this  equation  and  the  three  terms  in  the  parenthesis.  The  first  is  a  "mean" 
term,  the  second  is  a  "contribution"  due  to  the  variability  of  the  conductivity  field 
and  the  third  is  the  log  conductivity  and  head  gradient  cross-covariance.  Gelhar 
[1993]  evaluates  this  expression  for  isotropic  log  conductivity  fields  using  spectral 
methods.  The  result,  which  is  independent  of  the  shape  of  the  correlation  function, 
is  surprisingly  simple  [Gelhar,  1993] 

<fj>=-olJd  (2.19) 

where  d  is  the  dimension  of  the  domain.  This  negative  correlation  between  head  gradi- 
ent and  log  conductivity  has  a  profound  effect  on  the  flow.  The  effective  conductivity 
may  be  defined  as  that  value  when  multiplied  by  the  spatial  mean  gradient  results 
in  the  observed  specific  discharge  (flow  rate  divided  by  discharge  area).  The  a  pri- 
ori effective  conductivity  is  the  effective  conductivity  estimated  from  Equations  2.18 
and  2.19.  With  d  =  2,  we  have  for  two-dimensional  conductivity  fields,  the  a  priori 
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effective  conductivity  is  simply  Ke  =  Kg.  Notice  that  the  effective  conductivity  can 
be  significantly  less  than  the  arithmetic  mean  conductivity  Ka  =  Kgexp[afnk/2]. 

The  definition  of  the  effective  conductivity  is  quite  general,  and  a  general  effective 
conductivity  expression  is  denoted  Ke  =  Kgc  where  c  is  the  "correction"  due  to 
the  aforementioned  processes.  Under  this  definition,  there  is  always  some  effective 
conductivity  value,  regardless  of  nonstationarities  in  boundaries  or  flow  domains,  but 
its  magnitude  may  depend  upon  several  factors.  There  is  no  general  a  priori  effective 
conductivity  estimator.  The  estimator  given  here  is  limited  to  stationary  head  and 
conductivity  fields. 

2.5.2    Eulerian  Velocity  Field 

The  relationship  between  the  Eulerian  conductivity  field  and  its  associated  Eule- 
rian velocity  field  has  been  one  of  the  more  thoroughly  examined  aspects  of  subsurface 
hydrology.  For  the  example  problem,  closed-form  expressions  are  sought  for  the  Eu- 
lerian, or  spatial,  mean  velocity  U,  the  Eulerian  velocity  covariance  Cua[r],  and  its 
point  distribution,  given  that  the  underlying  conductivity  field  is  large,  stationary, 
lognormally-distributed,  exponentially  correlated,  and  characterized  by  parameters 
H\nk,  cfnfc'  and  Afnfc,  and  subject  to  the  conditions  prescribed  for  the  model  aquifer 
(see  Section  2.1  and  Figure  2.1). 

From  the  development  of  the  effective  conductivity  relationship  in  the  preceding 
section,  an  expression  for  the  mean  Eulerian  velocity  is 

KQcJ 

U  =  ~T-  (2-2°) 

where  c  =  1  for  a  two-dimensional  isotropic  aquifer.  The  c  is  retained  to  emphasize 
that  this  expression  is  not  the  "consistent  first-order  approximation"  of  the  velocity 
that  would  result  from  dropping  the  a?nk  and  <  fj  >  terms  in  Equation  (2.18). 
For  the  two-dimensional  isotropic  aquifer,  which  is  the  model  case,  these  two  terms 
happen  to  be  equal  in  magnitude  and  opposite  in  sign,  and  the  consistent  first-order 
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expression  and  effective  expression  are  the  same.  This  damping  effect  of  the  negative 
head  gradient  and  log  hydraulic  conductivity  is  largely  responsible  for  the  robustness 
of  these  perturbation  expansions.  The  truncated  terms  are  not  negligible  in  the  sense 
their  magnitude  is  small.  For  instance,  for  a2nk  =  1,  their  magnitude  is  one-half  that 
of  the  mean.  This  brings  us  to  an  important  philosophical  point.  A  primary  usefulness 
of  these  expansions  is  to  identify  the  dominant  "components"  of  the  process.  A  purely 
"consistent"  approach  could  lead  to  neglect  of  some  important  process.  For  instance, 
consider  a  strict,  first-order  expansion  of  Equation  (2.16).  Equation  (2.17)  would 
contain  no  f2  term,  and  Equation  (2.18)  would  contain  no  afnk  term.  The  resulting 
effective  expression  would  seriously  underestimate  the  specific  discharge,  since  the 
contribution  of  the  conductivity  variability  would  "carry  no  weight"  with  respect  to 
the  other  contributing  processes.   Analysis  of  the  neglected  higher  order  terms  is 
necessary  for  a  more  complete  understanding  of  the  process.  The  general  goal  of  the 
approach,  however,  is  to  identify  the  dominant  subprocesses  that  contribute  to  the 
larger  process,  and  the  relative  importance,  or  weight,  of  each  subprocess. 

Expressions  for  the  Eulerian  velocity  covariance  for  the  model  system  have  been 
given  by  various  researchers,  including  Graham  and  McLaughlin  [1989a]  and  Rubin 
[1990].  Shapiro  and  Cvetkovic  [1988]  develop  an  inverse  velocity  covariance  function 
using  a  perturbation  expansion  of  the  inverse  flow  velocity  and  spectral  methods. 
Consider  once  again  the  perturbed  Darcy's  law  expression  given  by  Equation  (2.16). 
Dividing  this  by  the  homogeneous  mobile  water  content,  neglecting  second-order 
terms  and  higher,  the  mean-removed  expression,  or  "velocity  perturbation,"  is 

u'  =  ^(f  +  j)  (2.21) 
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If  these  processes  are  stationary  stochastic  processes,  they  can  be  represented  repre- 
sented by  a  Fourier-Stieltjes  integral  (see  e.g.,  Bakr  et  al.  [1978]) 

poo 

u'=       exp[iM-x]dZu»[M]  (2.22) 

J  OO 

where  i  —  \/-T,  M  is  the  wave  number  vector,  and  dZu<  is  the  complex  random 
amplitude  of  the  velocity  perturbation.  The  Fourier-Stieltjes  representation  of  Equa- 
tion (2.21)  is 

dZu,  =  ^-(dZf  +  dZj)  (2.23) 

u 

Compare  this  expression  to  that  derived  by  Shapiro  and  Cvetkovic  [1988]  (their  Equa- 
tion (A4)) 

dZl/v,  =  -Yj{dZf  +  dZj)  (2.24) 

The  spectra  that  result  from  these  complex  random  amplitudes  are  structurally  the 
same,  differing  by  the  square  of  the  mean  values  by  which  they  are  multiplied  (e.g.,  by 
(KgJ/9)2  on  the  one  hand  and  (6/(KgJ))2  on  the  other).  In  general,  the  correlation 
structure  of  the  inverse  of  a  random  function  is  not  the  same  as  that  of  the  random 
function  itself.  The  first-order  series  approximation  of  the  inverse  velocity  is  linear  in 
the  velocity  itself  (e.g.,  l/v  ~  1  —  v ),  and  thus  the  velocity  correlation  characteristics 
are  recovered  as  an  approximation  of  the  inverse  velocity  correlation  characteristics. 
For  the  case  of  the  exponentially  correlated  isotropic  conductivity  field,  the  longitudi- 
nal velocity  covariance  of  a  two-dimensional  field  is  given  by  (Shapiro  and  Cvetkovic 
[1988]) 

Cu[s]  =  3alk(KgJ/6)2  (exp[-s](s-2  +  3*-3  +  3s~4)  +  s"2/2  -  3s~4)  (2.25) 
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where  s  =  x/\\nk-  The  corresponding  expression  for  a  three-dimensional  field  is 

Cu[s]  =  Salk(KgJ/6)2  (exp[-s](s-2  +  5s"3  +  12s"4  +  12s"5)  +  s"3  -  12s"5) 

(2.26) 

Equation  (2.25)  provides  a  superior  match  to  the  data  in  comparison  to  Equa- 
tion (2.13)  (see  Figure  2.17).  However,  the  expression  for  the  Eulerian  covariance 
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Figure  2.17:  Comparison  of  velocity  covariance  expression  Equation  (2.25)  to  the  ap- 
proximation Equation  (2.13)  for  ofnfc  —  1.  Covariances  are  presented  as  dimensionless 
correlation  functions. 

given  in  Equation  (2.13)  captures  the  dominant  features  of  the  velocity  covariance 
structure,  including  the  correlation  length  and  the  point  velocity  variance  of  Equa- 
tion (2.25)  in  a  simple  and  familiar  form. 

As  noted,  however,  Equation  (2.13)  exhibits  some  systematic  differences  from  the 
observed  behavior.  The  data  show  a  rapid  initial  decay,  followed  by  a  long-range 
persistence  (see  Figure  2.17).  This  behavior  is  examined  with  an  expression  adapted 
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from  Dagan  [1984] 

o  r  i       2    (r,    M     dR]nk4>[r]  d^r}\ 

Rua[r]  =  ^l„fc  (  R\nk[r\  ^  J  l^'J 

where  R\nk  is  the  log  conductivity  correlation  function,  R\nk<t>  is  the  log  conductivity 
and  head  cross-correlation  function,  and  74,  is  the  head  variogram.  From  this  ex- 
pression, glean  that  the  velocity  covariance  is  a  process  comprised  of  three  dominant 
subprocesses,  namely  R\nk,  R\nk<t>,  and  7^.  At  small  displacements,  all  three  of  these 
subcomponents  exhibit  rapid  trend  towards  zero  (see  Figure  2.18).  The  R\nk<t>  and  7^ 
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Figure  2.18:  Components  of  the  velocity  covariance  function. 

expressions  exhibit  effects  that  persist  well  beyond  that  of  the  log  conductivity  field. 
The  long-range  correlation  effects  of  the  head  field  in  the  mean  direction  of  flow  are 
cannot  be  precisely  predicted  with  the  simple  exponential  function  Equation  (2.13). 
In  order  to  recover  the  proper  correlation  length,  the  approximation  must  system- 
atically over-predict  the  velocity  correlation  at  small  separations.  The  "cost"  of  the 


47 

approximation,  however,  will  realize  significant  returns  in  the  form  of  simplicity,  as 
shown  in  subsequent  sections. 

Previous  works  have  focused  upon  the  moments  of  the  Eulerian  velocity  field  that 
result  from  mean  gradients  imposed  upon  heterogeneous  conductivity  fields.  These 
works  shall  be  extended  slightly  with  a  hypothesis  about  a  specific  case  that  allows 
a  variety  of  insights  to  be  drawn. 

Hypothesis:  The  point  velocity  distribution  of  an  Eulerian  velocity  field  induced 
by  an  stationary  mean  gradient  J  across  a  lognormally  distributed,  exponentially 
correlated  log  hydraulic  conductivity  field  characterized  by  parameters  n\nk,  <?\nk' 
and  Xfn  k  and  saturated  with  mobile  water  at  a  constant  and  homogeneous  content  6 
is  lognormal  with  approximate  parameters 


fJ-\nu  =  In 

and 


KgJC 


(2.28) 


alu  =  \n[l  +  balk]  (2.29) 

(see  Appendix  A  for  details  of  development) .  A  fundamental  assumption  in  addition 
to  the  lognormality  hypothesis  is  that  the  "real  space"  moments  of  u  are  \iu  =  U 
and  o\  =  U2bafDk.  Equation  (2.13)  is  adopted  as  an  approximation  of  the  velocity 
covarince  function. 

In  lieu  of  a  rigorous  proof,  the  results  of  a  numerical  experiment  in  which  statisti- 
cally independent  velocities  are  drawn  uniformly  in  space  from  simulated  flow  fields 
are  presented  (see  Figure  2.16).  The  normalization  of  the  velocities  in  Figure  2.16 
than  that  of  previous  figures  in  that  the  log  velocities  are  normalized  by  the  sample 
mean  and  standard  deviations  to  emphasize  the  lognormality  of  the  velocity  distribu- 
tion. Normalization  by  the  moments  predicted  by  Equations  (2.28)  and  (2.29)  would 
have  resulted  in  a  lognormal  distribution  function  that  deviates  from  the  standard 
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normal.  This  indicates  that  the  moments  nu  =  U  and  a\  —  U2bafnk  may  require  some 
further  investigation.  The  adequacy  of  these  estimates  are  discussed  in  subsequent 
sections. 

2.5.3    Streamtube  Weights 

The  reference  velocity  field  is  now  to  be  subdivided  into  streamtubes.  Each 
streamtube  may  be  considered  to  carry  some  "weight"  in  relation  to  other  streamtubes 
for  a  given  collection.  Consider  the  partitioning  of  the  flow  field  into  N  streamlines 
separated  by  an  equal  distance  along  the  IP.  This  is  an  equal-area  weighted  parti- 
tioning scheme,  and  results  in  N  sub-IPs  of  equal  area,  but  dissimilar  discharge.  The 
discharge  through  sub-IP  i  is 


where  Aa  =  A/N.  The  relative  flow- weight  of  this  streamtube  i  is  given  by 


An  equal-flow  partitioning  scheme  results  in  in  N  sub-IPs  of  equal  discharge,  but 
of  dissimilar  area.  The  discharge  through  sub-IP  i  is 


di  =  9AaUi 


(2.30) 


Wfi  =  Ui/U 


(2.31) 


d  =  9AciiVi 


(2.32) 


=  Q/N 


Solving  for  the  area  yields 


(2.33) 


The  relative  area-weight  of  this  streamtube  i  is 


Wai  =  U/Vi 


(2.34) 
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The  perceptive  reader  will  notice  a  slight  notational  change  between  the  two  parti- 
tioning schemes.  A  lowercase  u  represents  velocities  drawn  uniformly  in  space  and  a 
lowercase  v  represents  velocities  drawn  uniformly  in  flow.  The  distinction  is  impor- 
tant. It  is  through  these  sub-IPs  that  water  enters  the  streamtubes  and  a  description 
of  the  sub-IP  properties  is  essential  to  streamtube,  and  thus  flow  field,  descriptions. 
If  from  a  collection  of  N  sub-IPs  the  random  selection,  or  drawing,  of  one  sub-IP  is 
equally  likely  as  any  other,  the  weight  assigned  to  any  one  observation  from  M  samples 
would  be  1/M.  Assume  that  all  of  the  sub-IPs  are  sampled.  The  average  sub-IP  dis- 
charge is  the  same  for  the  equal-area  and  equal-flow  partitioning  schemes,  and  equal 
to  d  =  8AU/N.  However,  the  discharge  variance  is  zero  for  equal-flow  and  nonzero 
for  equal-area,  when  the  velocity  field  is  heterogeneous.  The  weights  expressed  in 
Equations  (2.34)  and  (2.31)  allow  the  properties  of  one  collection  of  streamtubes  to 
be  mapped  to  the  other  collection  of  streamtubes.  These  weights  combined  with  the 
assumed  values  of  the  Eulerian  velocity  field  permit  derivation  of  properties  of  the 
initial  velocities  of  the  equal-flow  partitioning  schemes.  The  expected  value  of  v  must 
equal  the  expected  value  flow-weighted  u.  Thus, 

<v>  =  V 

u 

=<  —u> 
U 

=  jjiA  +  l*u)  (2.35) 
=  U(l  +  balk) 

=  ^(1  +  W„fc) 
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The  variance  of  v  is  given  by  var[u]  =< v2  >  -  <v>2.  The  expected  value  of  v2  must 
equal  the  expected  value  of  flow- weighted  u2.  That  is 

2  ^2 

<vz>  =<  —u  l> 


=  ^exp[3/iinu  +  ^a2nu] 

=  -  exp[3(^,nu  +  a2nu/2)  +  Sa2n  J 


(2.36) 


=  U2(l  +  bclk) 
Thus,  the  variance  of  v  is  given  by 

™[v}  =  V2bo2nk  (2.37) 

It  is  worthy  to  note  that  these  results  do  not  follow  from  an  assumption  that  v  —  u2/U. 
That  is,  the  velocity  along  the  Lagrangian  trajectory  is  not  simply  the  flow-weighted 
velocity  along  an  Eulerian  trajectory. 

Hypothesis:  The  point  velocity  distribution  along  a  flow-weighted  Lagrangian 
trajectory  in  the  model  stationary  Eulerian  velocity  field  is  lognormal  with  approxi- 
mate parameters 


Mint;  =  ^ 

and 


KgJc 

e 


(2.38) 


a2nv  =  \n[l  +  ba2nk}  (2.39) 

Moreover,  the  process  is  assumed  second-order  stationary  in  displacements  along 

the  mean  direction  of  flow,  that  is,  at  control  planes,  and  is  characterized  by  the 
covariance  function 

Cvf[r]  =  V2ba2nk  exp[-6r/Alnjfc]  (2.40) 
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where  r  is  the  lag  along  the  mean  flow  direction.  Again  this  follows  from  the  require- 
ment that  the  distribution  match  the  assumed  "real  space"  v  moments  fiv  =  V  = 
U(l  +  ba?nk)  and  o2v  =  V2ba?nk.  The  correlation  structure  of  the  FW/LT  velocities  is 
similar  to  that  of  AW/ET  (see  Figure  2.12). 

Employing  these  velocity  distribution  hypotheses,  it  is  now  possible  to  derive  some 
relationships  between  the  different  flow  field  partitioning  schemes.  From  observation, 
it  is  inferred  that  the  velocities  along  area-weighted  Eulerian-trajectories  and  flow- 
weighted  Lagrangian-trajectories  are  second-order  stationary  in  displacements  along 
the  mean  flow  direction  (recall  Figures  2.10  and  2.11).  Having  adopted  or  postulated 
expressions  for  these  first  two  moments,  it  is  now  possible  to  derive  expressions  for 
the  nonstationary  velocity  means  and  variances  of  the  area- weighted  Lagrangian  and 
the  flow-weighted  Eulerian  streamtubes. 

2.5.4    Control  Plane  Oriented  Velocity  Moments 

Cvetkovic  et  al.  [1996]  postulated  a  semi-analytical  estimator  for  a  nonstation- 
ary mean  velocity  along  area-weighted  Lagrangian-trajectory  observed  at  control 
planes.  We  would  like  to  work  directly  from  the  assumed  stationary  flow- weighted 
Lagrangian-trajectory  velocity  field.  The  general  expression  for  the  expected  value 
of  the  area- weighted  Lagrangian-trajectory  velocity  is 


Resolve  the  stationary  mean  velocities  into  the  mean  of  the  flow-weighted  Lagrangian- 
trajectory  field  and  zero-mean  perturbation. 


(2.41) 


va[x]  =< 


U 


V{l  +  v'[x})> 


(2.42) 


V(l+v'[0]) 
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Expand  the  denominator  in  a  series  to  second  order 

va[x]  =<  U{1  -  v'[0]  +  v'[0]2){l  +  v'[x})  > 

=  U(l  +  &  "  ^)  (2-43) 

=  U(l  +  ba2nk(l-exp[-br/Xlnk})) 
In  a  similar  fashion,  an  expression  is  derived  for  the  nonstationary  mean  velocity 
along  a  flow-weighted  Eulerian-trajectory.  In  this  case,  it  is  the  stationary  Eulerian- 
trajectory  velocity  field  that  serves  as  a  basis.  The  mean  of  the  perturbation  expan- 
sion is  the  Eulerian  mean  velocity.  It  is  interesting  to  note  that  this  result  requires 
no  series  expansion. 

uf[x]  =<-^-u[x}> 

=  ±<U(l  +  u'[0])U(l  +  u'[x})> 

=  17(1  +  £=M) 

=  U(l  +  ba2nkexp[-br/\lnk]) 
The  nonstationary  velocity  variances  along  the  respective  cross-weighted  trajec- 
tories are  derived  in  a  similar  fashion.  The  general  expression  for  the  area- weighted 
Lagrangian-trajectory  velocity  variance  is 

var[ua[x]]  =<-^-(v[x]  -  va[x])2  >  (2.45) 
v[0\ 

This  equation  may  be  re-written 

v&r[va[x}]  =<  -^v[x]2  >  -va[x}2  (2.46) 
v[0\ 

Resolve  the  first  term  of  the  preceding  equation  into  perturbation  expressions 
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Expanding  the  denominator  in  a  series,  taking  the  expected  value,  and  neglecting 
third-order  and  higher  moments  yields 

<  ^-Axf  >=  UV(1  +  2ba2nk{\  +  exp]  -  bx/Xlak}))  (2.48) 
v[0\ 

In  order  to  improve  the  quality  of  the  estimate,  the  right-hand  side  (RHS)  of  the 
preceding  equation  is  required  to  match  the  known  extreme  values  of  the  left-hand 
(LHS)  side.  The  LHS  at  X  =  0  is 

<-T7^[0]2>  =  U  <v[0}> 
v[0]  1  J  1  J  (2.49) 

=  UV 

The  RHS  at  X  =  0  is 

RHS[0]  =  UV{\  +  2balk(l  -  exp[-W/Alnfc])) 

(2.50) 

=  UV 

The  statistical  independence  of  the  initial  velocity  and  the  velocities  as  x  goes  to 
infinity  are  exploited  to  calculate  the  LHS  as  x  — >  oo 

LHS[oo]  =<-^j><w[oo]2> 

=  ol  +  V2  (2-51) 
=  UV(l  +  bo?nk)2 

The  RHS  as  x  ->  oo  is 

RHS[oo]  =  UV(1  +  2bafnk(l  -  exp[-6oo/Alnfc])) 

=  1^(1  +  2^) 

The  RHS  is  deficient  a  term  {ba2nk)2.  Rewrite  the  RHS  thus 


(2.52) 


RHS  =  UV(1  +  (2ba2nk  +  (ba2nk)2)(l  -  exp[-6r/Alnfc]))  (2.53) 
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Substituting  Equations  (2.53)  and  (2.43)  into  Equation  (2.45)  yields 

vax[va[x]]  =  UV(l  +  (2balk  +  (kr2nfc)2)(l  -  exp[-6r/Alnfc])) 

(2.54) 

-(U(l  +  ba?nk(l-exp[-br/\lnk})))2 
Further  expansion  of  this  expression  offers  no  intuitively  obvious  simplifications  or 
insights  in  its  entirety.  However,  this  expression  may  be  approximated  by  dropping 
a  few  terms  thus 

varK[x]]  =  V  Vn*  "  U2balk{2balk  +  (6a2nfc)2)  exp[-6r/Alnfc]  (2.55) 

Two  unconventional  techniques  were  employed  in  this  derivation  that  shall  now  be 
defended.  The  first  is  requiring  the  series  expansion  to  match  the  known  endpoints. 
Simple  expressions  that  condense  the  system  behavior  into  an  compact  form  that 
allows  the  relative  contribution  of  the  component  processes  to  be  easily  grasped  are 
sought.  The  series  expansion  and  subsequent  truncation  is  itself  an  approximation. 
The  systematic  error  introduced  by  "missing"  a  known  endpoint  is  distracting.  In  the 
case  of  the  velocity  covariance,  the  shape  of  the  correlation  function  was  "sacrificed" 
for  simplicity,  but  the  correlation  length  was  retained.  Had  the  simplified  expression 
neither  matched  the  shape  of  the  data  nor  reproduced  the  correlation  length,  the 
expression  would  have  been  somewhat  more  questionable.  In  this  case,  the  expression 
is  required  meet  the  observed  endpoint  values  at  the  "cost"  of  the  rigor  of  the  series 
expansion  and  truncation. 

The  second  simplification  is  the  ad  hoc  dropping  of  terms  to  simplify  the  expres- 
sion. Again,  simple  expressions  that  allow  the  relative  contribution  of  the  different 
processes  be  easily  grasped  are  sought.  It  would  be  difficult  to  hypothesize  Equa- 
tion (2.55)  a  priori  as  an  empirical  expression.  A  primary  function  of  the  series 
expansion  is  to  indicate  the  dominant  processes  and  give  an  estimate  of  their  relative 
contributions.  This  hybrid  approach  is  useful  for  the  purposes  here.  These  techniques 
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are  not  a  panacea,  however,  and  the  approximations  have  only  been  tested  for  the 
log  conductivity  variances  tested  in  this  work. 

The  nonstationary  flow- weighted  Eulerian-trajectory  velocity  variance  may  be  de- 
rived in  a  similar  fashion.  These  results  are  presented  directly.  The  complete  expres- 
sion is 

var  [u,  [*]]  =  t/2  (l  +  ba2Xnk  +  [2bafnk  +  3  {ba2nk)2  +  (ftojU)3)  (1  -  exp  [-6r/Alnfc])) 

-(L/(l  +  6a2n,exp[-6r/Aln,]))2 

(2.56) 

and  an  approximation  of  the  above  obtained  by  neglecting  a  few  terms  is 

var  [uf  [x\]  =  U2bo?nk  (l  +  (2bo2nk  +  {balk)2)  exp  [-6r/Alnfc])  (2.57) 

2.5.5    Control  Time  Oriented  Velocity  Moments 

The  displacement-invariant  velocity  covariances  for  u  and  v  permitted  derivation 
of  expressions  for  the  first  two  moments  of  the  respective  velocities  for  the  four  dif- 
ferent trajectory  collections.  However,  only  the  area-weighted  Lagrangian-trajectory 
collection  appears  to  be  stationary  in  time.  This  is  perhaps  not  intuitively  obvious, 
as  the  velocity  is  clearly  not  stationary  in  displacements  along  the  mean  direction 
of  flow  (e.g.,  Equations  (2.43)  and  (2.54))  because  there  is  an  implicit  distribution 
of  times  associated  with  the  establishment  of  a  reference  injection  plane.  After  long 
travel  distances  and  times,  the  area-weighted  Lagrangian-trajectory  point  velocity 
statistics  are  equal  to  the  flow-weighted  Lagrangian-trajectory  point  statistics.  At 
the  injection  plane,  the  area-weighted  Lagrangian-trajectory  point  velocity  statistics 
are  identically  equal  to  the  area-weighted  Eulerian-trajectory  point  velocity  statistics. 
These  relationships  can  be  exploited  to  derive  equations  that  provide  a  few  insights 
about  how  properties  change  in  Lagrangian  time,  and  to  derive  an  expression  for  a 
Lagrangian-time  dependent  velocity  covariance  function. 
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The  area- weighted  velocity  at  the  IP  is  U  by  definition.  Since  the  area- weighted 
Lagrangian-trajectory  velocity  is  stationary  in  time,  its  mean  in  time  is  U  regardless 
of  position.  At  long  travel  times,  the  point  velocity  is  stationary  in  displacements  as 
well.  For  constant  means  in  space  and  time,  the  following  relationship  holds 

?/>-(*jr  iif  (2-58) 

This  expression  states  that  the  harmonic  mean  velocity  in  space  equals  the  arithmetic 
mean  velocity  in  time.  Returning  now  to  the  model  aquifer,  invoke  the  hypothesis  of 
lognormality  of  the  velocities  along  either  trajectory.  The  harmonic  mean  of  v[x]  is 
given  by 

Vh  =  exp[/iin„  -  (Tint,/ 2] 

=  exp[/iin„]  •  exp[<7in„]1/2 

  1  (2.59) 

=  U 

Similarly,  the  harmonic  mean  of  u[x]  is  given  by 

Uh  =  1^2  (2-6°) 

Although  derived  by  different  means,  this  result  is  equal  to  the  result  given  by  Shapiro 
and  Cvetkovic  [1988].  These  harmonic  means  are  the  asymptotic  mean  velocities  of 
the  two  trajectories  in  time.  The  initial  time-dependent  velocities  are  identical  to 
the  initial  space-dependent  velocities,  an  artifact  of  starting  all  of  the  trajectories 
at  the  injection  plane  at  an  implied  time  t  =  0.  This  indicates  that  the  velocity 
covariance  in  time  along  the  area- weighted  Lagrangian  trajectory  is  Cva[t  =  0]  = 
U2bo~\nk.  As  mentioned,  at  large  times,  the  statistical  properties  of  the  Lagrangian- 
trajectories  converge.  Making  a  change  of  variables  x  =  Vht  —  Ut'm  Equation  (2.40) 
and  substituting  the  known  Cva[t  =  0]  for  the  value  of  Cvf[r  =  0]  in  yields  an 
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expression  for  the  stationary  area-weighted  Lagrangian-trajectory  velocity  covariance 

Cva[t]  =  U2balkexp[-bUt/Xlnk]  (2.61) 

Similar  expressions  have  been  derived  before  by  other  researchers  (see  e.g.,  Dagan 
[1989]),  but  the  conceptual  development  is  quite  different.  Rather  than  seek  to  ap- 
proximate the  covariance  function,  an  approach  to  quantify  it  directly  from  properties 
of  the  system  is  demonstrated.  Moreover,  the  quantity  U  is  perceived  to  be  the  har- 
monic mean  of  the  Lagrangian-trajectory  velocity,  and  that  it  is  quantitatively  equal 
to  the  area-average  velocity  of  the  flow  field.  This  is  an  important  point. 

The  mean  velocity  and  velocity  variance  along  the  flow-weighted  Lagrangian- 
trajectory  is  found  in  a  manner  similar  to  those  in  space,  but  there  are  some  concep- 
tual issues  that  should  be  addressed.  It  is  the  collection  of  area-weighted  Lagrangian- 
trajectories  that  exhibits  a  stationary  velocity.  The  mean  velocity  is 

Vft[t]=<^v[t}>  (2.62) 

Expand  the  velocities  into  a  mean  and  a  zero-mean  perturbation.  The  mean  here  is 
the  time  average  mean,  or  the  harmonic  mean  of  the  space  average. 

«/<[<]  =  ^<U(l+v'[0})U(l  +  v'[t\)> 

=  U  (1  +  balkexp[-bUt/\lnk}) 
This  expression  requires  no  series  expansion.  For  the  control  plane  oriented  velocity 
statistics,  the  mapping  was  from  flow- weighted  to  area- weighted  trajectories.  Equa- 
tion 2.63  is  identical  in  form  to  Equation  2.44  with  the  substitution  r  =  Ut.  Recognize 
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the  quantitative  similarity  and  write  velocity  variance  equations  directly  The  com- 
plete expression  is 

var  [vf  [t]}  =  U2  (l  +  ba2nk  +  (26a,2n,  +  3  (ba2nk)2  +  {ba2nk)3)  (1  -  exp  [-bUt/ Xlnk})) 

-  (U  (l  +  ba2nkexV[-br/\lnk\))2 

(2.64) 

and  the  approximation  is 

var[v/[t]]  =  U2bafnk  (l  +  [2ba2nk  +  {ba2ak)2)  exp[-bUt/Xlnk})  (2.65) 

Calculation  of  the  mean  velocities  along  the  Eulerian  trajectories  is  problematic 
in  the  absence  of  a  stationary  velocity  covariance  function.  However,  note  that  the 
area-weighted  velocities  start  at  U  and  that  the  flow-weighted  velocities  start  at  V 
and  both  progress  to  U/(l  +  bafnk).  Recognizing  the  similarity  shared  by  the  different 
mean  time  estimators,  write  approximate  expressions  directly 

uat  [t]  =  ,Ji2    (1  +  ba2nk  exp  [-bUt/X{nk])  (2.66) 

and 

W  =  T^TT-  I1  +  (2^infc  +  *»lk)  (^P  [-Mtt/Ain*]))  (2-67) 

1  +  0<7ln  fc 

Again,  devoid  of  a  stationary  covariance  function  for  the  Eulerian  velocity  field  in 
time,  intuitive  arguments  lead  us  to  anticipate  an  asymptotic  u  variance  in  time  to 
be 

2        t  r2  l  2 
aut  =  Uhba\nk 

(     U     \\  2  (2-68) 

The  initial  velocity  variances  are  known  from  the  covariance  functions  in  space.  It  is 
possible  to  anticipate  forms  of  the  nonstationary  variance  equations  similar  to  those 
found  previously,  and  include  them  for  the  sake  of  symmetry  and  completeness.  The 


59 


area-weighted  Eulerian-trajectory  time-dependent  velocity  variance  approximation  is 
varK[i]]  =  (y^-J  Wnfc  (1  +  {^<k  +  {b<k))  exp  [-bUt/X^])  (2.69) 


The  flux-weighted  Eulerian-trajectory  time-dependent  velocity  variance  approxima- 
tion is 
var[uy[£]]  = 

balk  (l  +  (36aL  +  3  (6a,2nfc)2  +  (&a2nfc)3)  exp  [-6t/t/Alnfc]) 

(2.70) 


U  " 

1  + 


2.5.6    Log  Conductivity  Moments 

The  first  two  moments  of  the  log  conductivity  along  those  trajectories  for  which 
the  moments  are  stationary  at  control  planes  located  along  the  mean  direction  of 
flow  or  at  control  times  are  sought.  Working  from  a  stationary  Eulerian  reference 
field,  the  mean  and  variance  along  the  AW/ET  are  fi\nk  and  o?nk,  respectively.  Infer 
that  the  FW/LT  is  stationary  as  well,  analogous  to  the  velocity  field  results.  The 
flow-weighted  mean  log  conductivity  perturbation  is  given  by 
u[0]      _  Kgex?[f]J(l+j)9 

U  J  9KeJ         1  (2.71) 

=<(l  +  f  +  f2/2  +  ...)(l+j)f>/c 

Expanding,  retaining  second-order  terms,  and  taking  the  expected  value  yields 

<!f1/>=4,(l-l/^)/c  (2.72) 
which,  for  two-dimensional  fields,  is  crfak/2.  Similarly,  the  variance  is 

<^o]/2>  _  <!ffl/>2=  {<k  _         _  l/d)?)  /c  (2  73) 

which,  for  two-dimensional  fields,  is  cr2nfc(l  —  ofnfc/4) 
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Finally,  note  that  the  AW/LT  log  conductivity  appears  stationary  in  time  as  does 
the  velocity.  Moreover,  it  appears  to  have  values  equal  in  magnitude  to  its  AW/ET 
counterpart,  as  does  the  velocity. 

Unlike  the  velocity,  expressions  for  the  trajectory-based  cross-covariance  of  the 
velocity  and  log  hydraulic  conductivity  that  are  necessary  for  developing  the  non- 
stationary  log  hydraulic  conductivity  moments  were  not  developed.  This  is  left  for 
future  work. 

2.6    Evaluation  of  Simulation  and  Theoretical  Results 

2.6.1  Mean  Velocities  in  Space 

Figure  2.10  presents  the  control  plane  oriented  mean  velocities  for  the  four  different 
collections  of  streamtubes.  The  mean  velocities  are  normalized  by  U  =  KgJ/6  =  1/10 
and  this  value  was  used  in  the  mean  velocity  expressions  as  well.  That  is,  the  parame- 
ter values  are  the  input  parameters  to  the  simulations.  Other  input  parameter  values 
are  6  =  3/8  and  ofnfc  =  {1/2, 1}  depending  upon  the  set  of  simulations.  The  deviation 
of  the  mean  of  area- weighted  Eulerian-trajectory  data  from  1  are  an  indication  of  the 
quality  of  the  effective  conductivity  expression. 

The  area-weighted  Eulerian-trajectory  mean  velocity  appears  stationary,  in  ac- 
cordance with  established  theory,  and  its  value  appears  to  be  well  predicted  by  the 
effective  conductivity  relationship.  The  flux-weighted  Lagrangian-trajectory  mean 
velocity  appears  stationary,  supporting  the  stationarity  hypothesis.  Moreover,  the 
value  of  the  mean  velocity  appears  to  be  well  predicted  for  both  variances  of  log 
conductivity.  Additionally,  the  nonstationary  mean  velocity  estimators  appear  to 
accurately  reproduce  the  transition  between  the  asymptotic  velocities. 

2.6.2  Velocity  Variances  in  Space 

Figure  2.11  presents  the  control  plane  oriented  velocity  variances  for  the  four  dif- 
ferent collections  of  streamtubes.  The  variances  are  normalized  by  the  area-weighted 
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Eulerian-trajectory  velocity  variance  estimated  from  the  input  parameters  to  the  sim- 
ulations. That  is,  by 

^  =  (V/»)Vhi  (2-74) 

As  with  the  mean  velocity  estimators,  the  variance  estimators  employ  only  the  sim- 
ulation input  parameters. 

Once  again,  the  area- weighted  Eulerian-trajectory  variance  exhibited  its  pre- 
sumed behavior,  and  supports  the  well  established  fact  that  the  Eulerian  velocity 
field  is  fairly  well  understood  in  space.  The  hypothesized  flow-weighted  Lagrangian- 
trajectory  velocity  variance  is  somewhat  higher  that  the  data  for  both  afnk.  A  possible 
reason  is  the  negative  correlation  between  log  conductivity  and  head  gradient  is  not 
weighted  properly  in  the  collection  of  trajectories  skewed  towards  higher  log  conduc- 
tivity values. 

The  other  two  estimators  also  suffer  from  this  bias  in  the  asymptotic  Lagrangian- 
trajectory  velocity  variance,  although  the  transition  between  initial  and  final  values 
for  both  seems  to  be  well  characterized. 

2.6.3    Mean  Velocities  in  Time 

Figure  2.13  presents  the  time-dependent  mean  velocities  for  the  four  different  col- 
lections of  streamtubes.  Again,  the  mean  velocities  are  normalized  by  U  as  estimated 
from  the  simulation  input  parameters.  Only  the  area-weighted  Lagrangian-trajectory 
mean  velocity  exhibited  a  constant  mean  behavior  in  time.  Its  estimated  harmonic 
mean  is  U,  and  the  estimate  seems  good,  although  one  might  argue  that  there  is  a 
slight  over-predictive  bias. 

The  estimated  harmonic  mean  of  Eulerian  trajectories  appear  to  have  a  stronger 
under  predictive  bias.  The  transition  of  the  flow-weighted  Lagrangian-trajectory 
appears  to  be  captured  better  than  those  of  the  Eulerian  trajectories.  It  should  be 
noted  that  the  Lagrangian  trajectory  transition  was  based  upon  an  expression  for 
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a  stationary  velocity  field,  whereas  the  Eulerian  trajectories  are  ad  hoc  estimators, 
based  upon  a  transition  between  known  end  points. 

2.6.4  Velocity  Variances  in  Time 

Figure  2.14  presents  the  time-dependent  velocity  variances  for  the  four  differ- 
ent collections  of  streamtubes.  The  variances  are  normalized  by  the  area-weighted 
Lagrangian-trajectory  velocity  variance  estimated  from  the  input  parameters  to  the 
simulations.  The  area- weighted  Lagrangian- velocity  appears  to  be  stationary  in  the 
variance.  The  data  seem  to  lie  slightly  above  the  estimated  value. 

The  over-prediction  of  the  initial  flux- weighted  velocity  variance  has  far  less  impact 
upon  the  time-dependent  velocity  variances.  This  combined  with  what  appears  to 
be  a  good  estimate  of  the  long-time  Eulerian-trajectory  velocity  variance  yields  a 
better  match  to  the  simulations  than  one  might  expect  from  the  means.  However, 
there  is  a  strong  initial  decline  in  the  variance  not  captured  by  the  simple  exponential 
expressions.  The  area-weighted  Lagrangian-trajectory  covariance  function  exhibits  a 
similar  rapid  decay.  The  early  time  behavior  is  likely  dominated  by  the  conductivity 
covariance.  As  the  trajectory  experiences  a  wide  range  of  conductivity  variability, 
this  effect  diminishes,  and  the  longer  range  correlations  related  induced  by  the  head 
field  begin  to  dominate. 

2.6.5  Log  Conductivity  Statistics 

The  predictors  of  the  stationary  log  conductivity  statistics  worked  well  for  both 
the  space-dependent  and  time-dependent  means  although  there  seemed  to  be  a  slight 
over-predictive  bias  for  the  flow- weighted  trajectory  (see  Figures  2.6  and2.8).  The 
AW/ET  log  conductivity  variance  was  "well  predicted"  in  that  the  turning-bands 
method  performed  well,  as  discussed  in  Sections  2.4.1  and  2.4.7.  The  log  conductivity 
variance  of  the  FW/LT  appeared  to  be  lower  than  that  of  the  AW/ET,  as  predicted  by 
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Equation  (2.73),  but  the  magnitude  of  under  prediction  was  larger  than  that  observed 
(see  Table  2.4). 

Table  2.4:  Comparison  of  observed  log  hydraulic  conductivity  variances  to  estimators. 
Observed  value  is  zero-slope  line  regressed  through  data.  The  decimal  representation 
of  the  fractional  value  is  included  for  convenience. 

<k  =  1/2 

parameter  predicted  observed 

AW/ETai2nJt[x]  1/2  0497" 

FW/LTcrj^x]  7/16=0.4375  0.484 

AW/LTajUM  1/2  0.498 

parameter  predicted  observed 
M/ETalk[x]            1  L002~ 
FW/LTafnk[x]          3/4  0.944 
AW/LTagJi]            1  1.000 

2.7    General  Conclusions 

In  an  steady  irrotational  flow  field,  a  fluid  parcel  trajectory  follows  a  hydrodynamic 
streamline.  The  flow  field  properties  as  observed  along  streamlines  can  differ  radically 
than  those  observed  "uniformly  in  space."  For  the  model  system,  a  stationary  Eu- 
lerian  log  conductivity  field,  this  work  identified  the  trajectory  collections  for  which 
certain  properties  appear  statistically  stationary  to  second  order.  For  displacements 
along  the  mean  flow  direction,  point  velocity  and  log  conductivity  appear  stationary 
along  streamlines  separated  by  equal  discharges.  For  displacements  in  travel  time, 
these  same  properties  appear  stationary  along  streamlines  that  are  selected  uniformly 
in  space.  The  stationarity  of  these  properties  facilitates  the  development  of  simple 
expressions  for  the  statistics  of  displacement  and  travel  times,  as  we  shall  show  in 
following  chapters. 

An  immediate  result  of  the  preceding  analysis  is  that  the  bulk  of  the  flow  field  in 
the  sense  of  total  discharge  travels  at  a  velocity  that  is  locally  greater  than  the  spatial 
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average,  and  this  disparity  increases  with  heterogeneity  of  the  underlying  conductivity 
field.  This  might  be  described  as  "preferential  flow."  A  practical  implication  for  the 
transport  of  kinetically  sorbing  or  interacting  solutes  is  that  the  local  velocity  for  the 
majority  of  the  flow  may  occur  in  a  regime  significantly  higher  than  that  estimated 
as  the  flow  field  average.  An  understanding  of  the  medium  heterogeneity  is  essential 
if  laboratory  experiments  are  to  be  "scaled  up." 

It  is  intuitively  obvious,  and  has  certainly  been  noted  before,  that  the  bulk  of  the 
flow  passes  through  an  increasingly  small  portion  of  the  total  "swept  volume"  as  the 
medium  heterogeneity  increases.  Again,  the  practical  implications  are  likewise  obvi- 
ous. Consider  some  aquifer  remediation  effort  that  involves  flushing  a  "contaminated 
zone"  with  some  fluid.  This  might  be  as  simple  as  the  upgradient  water  in  the  case  of 
pump-and-treat  or  as  extravagant  as  an  surfactant-enhanced  microemulsion.  In  het- 
erogeneous flow  systems,  the  bulk  of  the  flushing  fluid  contacts  a  disproportionately 
small  volume  of  the  total  swept  volume.  For  contaminants  that  tend  to  reside  in  areas 
of  high  conductivity,  the  serves  to  enhance  the  efficacy  of  remediation  efforts.  Con- 
taminants that  tend  to  reside  in  areas  of  low  permeability,  conversely,  would  be  more 
difficult  to  remove.  While  these  conclusions  may  be  drawn  from  "common  sense," 
application  of  the  concepts  developed  here  may  help  to  better  understand  these  issues 
quantitatively. 


CHAPTER  3 
NONREACTIVE  SOLUTE  TRANSPORT 

3.1  Introduction 

The  previous  chapter  saw  the  derivation  of  the  relationship  between  the  velocities 
that  lie  along  Lagrangian  trajectories  and  the  parameters  that  describe  a  simple 
model  Eulerian  flow  field.  These  Lagrangian  relationships  are  employed  to  derive 
a  few  functions  that  characterize  global  measures  of  solute  transport:  namely  the 
statistics  of  mass  displacement  in  time  and  the  statistics  of  mass  arrival  times  at 
control  planes.  The  frameworks  are  are  that  of  Dagan  [1982b]  and  that  of  Shapiro 
and  Cvetkovic  [1988],  respectively.  The  recovery  of  a  few  well-known  results  are 
anticipated,  in  addition  to  some  that  are  new.  However,  the  approach  is  novel  in  that 
the  work  is  based  upon  trajectory-based  statistical  properties  of  the  flow  field. 

3.2    Injection  Mode  and  Streamtube  Selection 

Demmy  et  al.  [1999]  noted  the  relationship  between  boundary  conditions  and 
sampling  of  streamlines.  A  clarification  of  this  relationship  and  a  generalization  the 
results  of  that  work  are  sought.  Unlike  the  previous  chapter,  the  Eulerian  trajecto- 
ries are  not  considered  in  the  context  of  an  actual  transport  trajectory,  as  they  are 
of  little  practical  importance.  However,  the  results  of  Shapiro  and  Cvetkovic  [1988] 
correspond  to  area-weighted  Eulerian-trajectories,  as  will  be  shown.  Therefore,  tra- 
jectory, Lagrangian  trajectory,  streamtube,  and  streamline  are  used  interchangeably. 

3.2.1    Uniform  Resident  Injection 

Consider  the  application  of  mass  at  a  constant  density  c0  across  the  entire  face  of 
the  injection  plane  (IP)  of  the  model  aquifer  (see  Figure  2.1).  Let  the  mass  occupy 
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some  small  Ax  that  is  invariant  with  displacements  along  the  IP.  Following  the  ter- 
minology of  Kreft  and  Zuber  [1978],  this  injection  mode  is  defined  to  be  a  uniform 
resident  injection.  Uniform  refers  to  the  homogeneous  mass  density  in  space  and 
resident  refers  to  the  volume  of  resident  fluid  into  which  the  solute  is  introduced. 
Subsequent  references  to  a  resident  injection  will  refer  to  a  uniform  resident  injection 
unless  otherwise  noted.  Now,  this  injection  mode  has  nothing  to  do  with  the  flow  field 
per  se.  That  is,  the  mass  "magically  appears"  and  occupies  a  fixed  volume  at  a  fixed 
concentration  without  being  influenced  by  local  water  fluxes.  Consider  the  flow  field 
to  be  divided  into  some  number  of  streamtubes,  say  N.  Moreover,  let  these  stream- 
tubes  carry  an  equal  area- weight  at  the  IP.  For  a  small  depth  of  injection  Ax,  this 
injection  mode  has  assigned  an  equal  mass  weight  to  each  streamtube  with  respect 
to  the  total  solute  mass.  The  total  mass  of  solute  is 

Mr  =  cQ6hLAx  (3.1) 

The  mass  weight  for  any  one  streamtube  i  is 

^  =  1  (3.2) 
M     N  v  ' 

This  value  is  identically  equal  to  the  area- weight  of  the  streamtube,  and  constant  for 
all  in  the  collection. 

Consider  now  dividing  the  streamtube  into  N  streamtubes  that  all  carry  the  same 
flow  weight.  The  mass  in  some  streamtube  i  is  given  by 

¥  =  b  (3.3) 
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3.2.2    Injection  in  Flux 

Consider  now  maintaining  the  IP  at  c0  for  some  brief  interval  At.  A  mass 

Mj  =  c08hLUAt  (3.4) 

enters  the  domain  and  the  different  streamtube  collections.  This  injection  mode  is  a 
uniform  injection  in  flux.  Uniform,  again,  refers  to  the  uniform  mass  density  and  "in 
flux"  refers  to  the  influent  water  that  carries  the  solute  into  the  flow  domain.  A  mass 

rrii  =  co0h^ViAt  (3.5) 

enters  equal-area  streamtube  i.  This  corresponds  to  a  mass  weight  of 


A  mass 


rrii  =  coOhciiViAt 


(3.6) 


Q  (3-7) 
=  c09h^At 

enters  equal-flow  streamtube  i.  This  corresponds  to  a  mass  weight  of  1/N. 

3.3    Integrated  Streamtube  Properties 

Let  us  return  to  the  notion  of  the  streamtubes  as  pure  hydrodynamic  entities.  It 
is  possible  to  discuss  properties  of  these  streamtubes  with  no  mention  of  injection 
modes.  Two  properties  hold  some  hydrodynamic  interest.  The  first  is  the  velocity 
integrated  over  time.  For  for  some  streamtube  i  this  is 


Xi[t]=  [  Vi[t]dt  (3.8) 
Jo 

A  velocity  integrated  over  time  is  a  "displacement"  with  a  dimension  of  length.  The 
second  the  inverse  velocity  integrated  over  some  distance.  For  some  streamtube  i  this 
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is 


(3.9) 


This  integral  has  a  dimension  of  time.  These  integrals  are  given  in  terms  of  the  x 
velocity  component  as  a  function  of  longitudinal  displacement  or  time.  It  is  suggested 
that  future  works  examine  the  details  of  a  more  general  approach,  e.g.,  displacements 
along  the  trajectory. 

3.3.1  Displacements 

Consider  a  model  flow  field  such  as  the  one  given  in  the  previous  chapter.  Quan- 
tification of  the  displacement  and  travel  time  statistics  for  area-weighted  and  flow- 
weighted  trajectories  is  sought.  Previous  results  indicate  that  the  velocity  for  a 
collection  of  area-weighted  streamtubes  is  stationary  in  time,  characterized  by  mean 
U  and  approximate  covariance  U2bafnk  exp[— bUt/X\nk\-  The  expected  value  of  Equa- 
tion (3.8)  is 


Although  the  velocity  of  the  flow-weighted  streamtube  collection  is  not  stationary  in 
time,  the  results  of  the  previous  chapter  are  employed  to  map  the  properties  of  the 
area-weighted  collection  to  the  flow-weighted  collection. 


Recognize  that  this  equation  may  be  re-arranged  to  be  the  integral  of  Equation  (2.63) 
in  time.   Substitute  the  RHS  of  Equation  (2.63)  into  the  preceding  equation  and 


<xa[i}>  =<  /  v[t]dt> 


(3.10) 


=  Ut 


(3.11) 
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integrate  to  yield 


<xf\i}>=  Ut  + X]nkalk{l-exp[-bUt/Xlnk}) 


(3.12) 


The  displacement  variance  for  the  area-weighted  streamtube  collection  is  given  by 
var[xa[t]]  =<xa[t]2>  -  <xa[t}>2 


3.3.2    Travel  Times 

Expressions  for  the  travel  time  statistics  are  sought,  in  analogy  to  those  sought  for 
displacement.  The  reference  collection  of  velocities  from  which  these  expressions  will 
be  derived  is  the  stationary  velocities  of  the  equal  flow  streamtubes,  or  the  equal-flow 
weight  Lagrangian  trajectories.  The  mean  travel  time  for  these  streamtubes  is  given 


(3.13) 


2  f \t  -  s)Clat[s}ds  -  (Ut)2 


by 


(3.14) 


x 


U 


The  mean  travel  time  for  the  area-weighted  streamtubes  is 


(3.15) 


=  -(x  +  Xlnka2nk{l  -  exp[-bx/Xlnk})) 


70 

The  arrival  time  variance  for  the  flow-weighted  collection  of  streamtubes  is 


Jo  v[x'\  Jo  v[x"\ 


Approximate  the  inverse  velocity  covariance  contained  in  the  first  term  of  the  RHS 
of  Equation  3.16  with  the  a  series  expansion  of  the  perturbation  equation 


1  1 

>  =<T7^7Z  ,  r    .,w.  7T-. 7TV> 


u[x']u[a;"]  V2{1  +  v'[x']){l  +  v'[x"]) 

=<  4(1  -  v'[x']  +  v'[xf  -...)(!-  v'[x"\  +  v'lx'f  -...)> 

r 

1 


V2  (3.17) 


=  ^(l  +  2a2v  +  Cvf[x',x"}+UOT) 


1   (,  ,  Cvf[x',x"Y 

U2  yL  +  y2 


where  HOT  is  the  collection  of  higher  order  terms.  Inserting  this  expression  into 
Equation  3.16  yields  an  expression  for  the  travel  time  variance  for  flux- weighted 
streamlines 

var[t/N]  =  ~f  (exp[-WA2nfc]  "  1)  j  (3-18) 

3.4    Solute  Transport 

Expressions  for  the  statistics  of  the  integrated  trajectory  properties  of  displace- 
ment and  travel  time  have  been  derived  in  the  absence  of  an  explicit  reference  to 
solute  transport  to  demonstrate  that  these  properties  exist  simultaneously  in  the 
same  flow  field.  However,  they  are  intimately  related  to  solute  transport.  Consider 
the  transport  of  an  infinitesimal  particle  whose  movement  is  constrained  to  follow  a 
streamline.  Moreover,  at  any  time  or  position  on  the  streamline,  the  particle  velocity 
is  identically  equal  to  that  of  the  underlying  velocity  field.  Its  displacement  from 
x  =  0  in  t  =  T  is  given  by  Equation  (3.8)  and  the  time  required  to  travel  from  x  =  0 
to  x  =  X  is  given  by  Equation  (3.9).  For  the  same  flow  field,  two  collections  of 
streamtubes  have  been  discussed,  namely  area-weighted  and  flow-weighted.  Placing 
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one  particle  into  each  area-weighted  streamtube  at  x  =  0  and  t  =  0  is  analogous 
to  the  a  uniform  resident  mass  injection.  Equations  (3.10)  and  (3.13)  express  the 
mean  particle  displacement  and  particle  displacement  variance  as  a  function  of  time. 
Equation  (3.15)  expresses  the  mean  particle  travel  time  at  different  control  planes. 

Similarly,  placing  one  particle  into  each  flow-weighted  streamtube  at  x  =  0  and 
t  =  0  is  analogous  to  the  a  uniform  injection  in  flux.  Equation  (3.12)  expresses 
the  mean  particle  displacement  as  a  function  of  time.  Equations  (3.14)  and  (3.18) 
express  the  mean  particle  travel  time  and  particle  travel  time  variance  at  different 
control  planes. 

In  addition  to  these  two  injection  modes,  there  are  hybrid  injection  modes.  For 
example,  consider  if  the  introduction  were  somewhat  random,  or  particle  mass  were 
some  deterministic  function  of  the  injection  velocity,  or  some  combination.  In  fact, 
it  is  left  to  the  imaginative  reader  to  conceive  of  different  injection  mode  possibilities 
and  the  physical  systems  to  which  these  correspond. 

3.5    Evaluation  of  Expressions 

In  the  previous  chapter,  a  numerical  experiment  used  to  test  certain  hypothesis 
about  the  relationship  between  Lagrangian  and  Eulerian  descriptions  of  flow  fields 
was  described.  Those  data  are  employed  here.  In  the  comparison  of  the  derived 
expressions  to  the  data,  the  functions  are  based  upon  the  simulation  input  parameters, 
and  all  normalization  is  carried  out  using  these  same  parameters.  Therefore,  the 
expressions  represent  a  priori  estimates  of  the  simulation  results,  based,  of  course, 
on  known  input  parameters. 

3.5.1  Displacements 

Figure  3.1  presents  the  model  estimates  of  mean  mass  displacements  for  the  two 
injection  modes.  The  flux  injection  preferentially  selects,  or  mass  weights,  high  veloc- 
ity areas  of  the  flow  domain,  and  the  propagation  of  the  center  of  mass  is  characterized 
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by  a  higher  velocity.  The  time- varying  mean  velocity  expression  for  the  flow- weighted 
Lagrangian-trajectories  (Equation  (2.63))  is,  in  fact,  an  expression  for  the  mean  ve- 
locity of  a  plume  resulting  from  a  uniform  flux  injection.  The  long-time  asymptotic 
slope  of  the  flux-weighted  mean  displacement  is  U,  as  is  easily  verified  by  taking  limit 
of  Equation  (2.63)  as  t  ->  oo.  The  mean  displacement  functions  predict  the  observed 
behavior,  and  also  provide  an  excellent  fit  to  the  data. 

The  displacement  variance  data  are  interesting  (see  Figure  3.2).  Injection  mode 
seems  to  have  little  effect  upon  the  displacement  variance  in  time.  This  is  not  to 
say  that  injection  mode  has  little  effect  upon  the  spatial  distribution  of  the  plume 
in  time.  At  a  given  time  t,  a  plume  injected  in  flux  will  have  traveled  a  greater 
distance,  on  average.  After  traveling  the  same  distance,  the  displacement  variance 
of  the  flux-injected  plume  will  be  less  than  that  of  the  resident  injection.  Thus, 
the  "spreading  as  function  of  mean  displacement"  is  less  in  the  case  of  the  flux- 
injected  solute.  The  area-weighted  trajectory  displacement  variance  over-predicts  the 
observed  displacement  variances,  and  the  long-time  slopes  appear  to  have  different 
values.  After  10  characteristic  times,  the  exponential  term  in  Equation  (3.13)  retains 
a  little  more  that  two  percent  of  its  original  value.  This  difference  does  not  seem  large 
enough  to  account  for  the  observed  discrepancy.  The  variance  of  the  stationary  area- 
weighted  Lagrangian-trajectory  velocity  field  in  time  appears  to  be  well  predicted. 
Thus,  the  problem  may  lie  with  the  rather  poor  description  of  the  short-time  and  long- 
time correlation  exhibited  by  the  observed  velocity  covariance  given  by  the  assumed 
model. 

3.5.2    Travel  Times 

Figure  3.3  shows  good  agreement  between  the  estimated  and  observed  mean  travel 
times.  The  injection  in-flux  shows  earlier  arrival  than  that  of  the  uniform  resident 
injection.  Uniform  resident  injection  weights  all  areas  of  the  IP  equally.  Large  areas 
that  contribute  little  to  the  overall  flow  receive  and  equal  mass  weighting  as  those 
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Figure  3.1:  Comparison  of  mean  displacements  for  the  two  injection  modes  or  stream- 
tube  collections.  The  flow-weighted  Lagrangian  trajectory  (FW/LT)  statistics  cor- 
respond to  a  flux  injection,  and  the  area- weighed  Lagrangian  trajectory  (AW/LT) 
statistics  correspond  to  a  uniform  resident  injection.  The  top  figure  contains  the 
aink  =  1/2  data  and  the  bottom  figure  contains  the  afnk  =  1  data. 
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Figure  3.2:  Comparison  of  displacement  variances  for  the  two  injection  modes  or 
streamtube  collections.  The  flow-weighted  Lagrangian  trajectory  (FW/LT)  statistics 
correspond  to  a  flux  injection,  and  the  area- weighed  Lagrangian  trajectory  (AW/LT) 
statistics  correspond  to  a  uniform  resident  injection.  The  top  figure  contains  the 
°fnfc  =  V2  data  and  the  bottom  figure  contains  the  ofnA.  =  1  data.  Variances  are 
normalized  by  Xikafnk. 
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areas  that  contribute  more  significantly  to  the  overall  flow.  However,  as  the  solute 
moves  through  the  domain  and  experiences  the  full  range  of  velocity  variability  along 
its  trajectory,  the  "slow  start"  is  damped  out  and  the  instantaneous  average  point 
velocity  for  all  for  the  particles  is  V  and  the  harmonic  mean  Vh  =  U.  The  travel  time 
variances  appear  to  be  better  matched  than  the  displacement  variances  (see  3.4).  This 
is  reasonable  in  light  of  the  better  apparent  match  of  the  assumed  velocity  covariance 
function  in  space  than  that  in  time.  The  uniform  injection  appears  to  exhibit  less 
non-linear  behavior  than  the  injection  in  flux.  The  uniform  injection  variance  is 
reasonably  expected  to  be  larger  than  that  of  the  flux  injection,  due  to  the  skewed 
mass-weighting  of  low-velocity  areas  by  the  uniform  resident  injection. 

3.6    General  Conclusions 

It  is  clear  that  the  injection  mode  effects  are  most  prevalent  in  the  "near  field" 
where  the  initial  velocities  and  travel  times  or  displacements  are  well  correlated.  Af- 
ter the  plume  travels  several  correlation  scales,  these  injection  mode  effects  diminish. 
The  near  field  and  the  near-to-far  transitional  scales  are  quite  interesting,  as  these 
intermediate  scales  are  on  the  order  of  sites  associated  with  certain  agricultural  set- 
tings such  as  some  feed  lots  and  waste  holding  facilities,  urban  land  uses  such  as  gas 
stations  and  dry  cleaners,  and  their  associated  remediation  efforts,  and  a  plethora  of 
"field  scale"  experiments.  These  scales  are  examined  in  the  perspective  of  the  "error" 
in  travel  time  estimate  associated  with  assuming  a  resident  injection  when  the  "real" 
injection  is  in  flux 

e[x]  =  Xinkafnk(l  -  exp[-bx/\lnk])/x  (3.19) 

given  by  subtracting  Equation  (3.15)  by  Equation  (3.14)  and  dividing  the  difference 
by  Equation  (3.14)  (see  Figure  3.5)  The  "near  field"  might  be  characterized  as  the  first 
five  to  ten  correlation  lengths,  where  the  injection  mode  effects  are  quite  prevalent, 
even  for  modest  log  conductivity  variability.  These  effects  diminish  less  rapidly  over 
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Figure  3.3:  Comparison  of  mean  arrival  times  for  the  two  injection  modes  or  stream- 
tube  collections.  The  flow- weighted  Lagrangian  trajectory  (FW/LT)  statistics  cor- 
respond to  a  flux  injection,  and  the  area-weighed  Lagrangian  trajectory  (AW/LT) 
statistics  correspond  to  a  uniform  resident  injection.  The  top  figure  contains  the 
a\nk  =  1/2  data  and  the  bottom  figure  contains  the  o^nk  =  1  data. 
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Figure  3.4:  Comparison  of  arrival  time  variances  for  the  two  injection  modes  or 
streamtube  collections.  The  flow-weighted  Lagrangian  trajectory  (FW/LT)  statistics 
correspond  to  a  flux  injection,  and  the  area-weighed  Lagrangian  trajectory  (AW/LT) 
statistics  correspond  to  a  uniform  resident  injection.  The  top  figure  contains  the 
aink  =  1/2  data  and  the  bottom  figure  contains  the  <rfnk  =  1  data.  Variances  are 
normalized  by  (t//(Ain*crinA:))2. 
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Figure  3.5:  Travel  time  estimate  error  as  function  of  displacement  associated  with 
assuming  a  uniform  resident  injection  when  injection  is  in  flux  for  0fnfc  =  {1/2, 1}. 

a  "transition"  region  between  10  and  30  correlation  lengths.  Beyond  30  correlation 
lengths,  the  difference  between  the  two  modes  is  less  than  a  few  percent. 

Much  of  the  Lagrangian  transport  literature  is  predicated  upon  a  uniform  resident 
injection,  but  the  so-called  consistent  first-order  travel  time  estimate  that  is  regularly 
employed  is  quantitatively  that  for  a  flux  injection.  For  large  displacements,  this 
error  is  small,  and  the  simplification  afforded  by  this  approximation  is  considerable 
(compare  Equation  (3.14)  and  Equation  (3.15)).  For  many  practical  applications, 
an  error  of  ten,  or  even  more,  percent  might  be  considered  negligible.  Therefore, 
it  is  not  suggested  that  this  approximation  be  abandoned,  but  rather  phrasing  such 
as  "inconsistent"  first-order  approximation  be  applied  in  future  works,  since  this 
Lagrangian  expression  is  inconsistently  based  upon  the  Eulerian  mean  velocity. 


CHAPTER  4 
REACTIVE  SOLUTE  TRANSPORT 

4.1  Introduction 

The  concepts  developed  in  the  previous  chapters  are  employed  to  analyze  reac- 
tive solute  transport.  The  focus  is  upon  global  measures  of  transport,  rather  than 
predictions  of  local  quantities.  Particular  processes  of  interest  are  injection  mode 
effects,  the  effects  of  correlation  of  the  local  "reaction  parameter"  to  the  log  con- 
ductivity field,  and  the  interplay  of  these  injection  mode  and  correlation  effects.  A 
suite  of  numerical  experiments  illustrates  these  effects.  Additionally,  expressions  are 
developed  for  temporal  moments  of  mass  breakthrough  curves  for  a  solute  subject  to 
linear  equilibrium  sorption  for  the  different  injection  modes  that  reflect  the  observed 
behavior  in  relatively  simple  closed-forms. 

Spatial  moments  of  sorbing  solutes  are  not  directly  treated  in  this  work.  The 
experimental  design  was  targeted  directly  towards  a  temporal  moment  analysis,  and 
this  design  decision  has  some  implications  for  spatial  moment  analysis.  These  decision 
and  implications  are  discussed,  in  hopes  that  future  works  might  be  guided  by  this 
analysis.  Some  of  the  properties  of  the  travel  time  dependent  reaction  parameter  are 
discussed,  as  are  some  of  the  implications  for  solute  transport. 

4.2    Temporal  Moments 

The  concepts  developed  in  the  preceding  chapters  are  applied  to  examine  the  effect 
of  injection  mode  upon  the  temporal  moments  of  the  breakthrough  curve  of  a  solute 
subject  to  linear  equilibrium  sorption.  A  uniform  mean  gradient  induces  a  steady 
irrotational  flow  in  a  uniformly  saturated  porous  medium  characterized  by  an  expo- 
nentially correlated  lognormal  conductivity  distribution.  This  flow  field  is  divided 
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into  collections  of  noninteracting  streamtubes.  Dispersion  at  the  streamtube  scale 
is  considered  to  be  negligible  in  its  effect  on  the  aggregate-scale  temporal  moments 
(e.g.,  Dagan  [1989]).  The  sorption  process  is  related  to  a  generic  reaction  parameter 
that  is  heterogeneous  and  may  be  correlated  to  the  log  conductivity  field. 

The  previous  chapters  demonstrated  the  relationship  between  streamtube  selec- 
tion strategies  (equal-area  or  equal-flow)  and  injection  mode  (uniform  resident  and 
in-flux).  The  reference  collection  is  equal-flow  streamtubes,  since  the  associated  ve- 
locity and  log  conductivity  statistics  are  stationary  in  displacements  along  the  mean 
flow  direction. 

The  primary  objective  of  this  chapter  is  to  examine  the  effect  of  injection  mode 
upon  the  temporal  moments  of  the  mass  breakthrough  curve  of  a  solute  subject  to 
linear  equilibrium  sorption  in  a  heterogeneous  flow  and  sorption  field.  A  brief  review 
the  conceptual  model  of  solute  transport  is  given.  General  expressions  are  presented 
for  the  first  two  temporal  moments  of  the  breakthrough  curves  associated  with  the 
uniform  injection  of  a  solute  both  in  the  influent  fluid  flux  and  in  the  resident  fluid. 
The  results  of  the  reactive  solute  transport  simulations  are  presented,  followed  by  the 
development  of  some  analytical  expressions  that  help  explain  the  observed  behavior. 

4.3    Conceptual  Overview  of  Solute  Transport 

Consider  a  general  multi-dimensional  flow  field  that  is  resolved  into  a  collection 
of  noninteracting  streamtubes.  Transport  along  an  individual  streamtube  is  charac- 
terized with  a  transfer  function  7[£;  x]  (see  Jury  and  Roth  [1990]  and  Cvetkovic  et  al. 
[1998]).  For  a  conservative  and  nonreactive  solute,  the  transfer  function  is  simply 

j[t;x]  =  S[t-T[x}]  (4.1) 

where  8  is  the  Dirac  delta  function  and  r  is  the  travel  time.  The  physical  interpreta- 
tion of  this  transfer  function  is  this:  a  particle  located  at  x  =  0  "arrives"  at  x  —  X 
at  time  t  =  r[X}.  We  conceive  of  this  travel  time  as  not  only  the  time  required  to 
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travel  to  some  point  x,  but  as  an  intrinsic  integrated  streamtube  property,  namely 
the  inverse  velocity  integrated  along  the  trajectory  (e.g.,  Equation  (3.9)). 

Recall  reaction  flow  path  concept  of  Cvetkovic  et  al.  [1998].  The  reaction  flow 
path,  parameterized  in  space,  is  defined  as 

where  P[x]  is  some  reaction  parameter  of  interest.  For  linear  equilibrium  sorption, 
P  is  a  linear  partitioning  coefficient.  The  reaction  flow  path  is  also  an  integrated 
streamtube  property.  As  such,  any  reference  to  /x  properties  in  general  are  broadly 
applicable  to  any  distributed  reaction  parameter,  perhaps  a  first-order  decay  coef- 
ficient or  a  Langmuir  sorption  parameter.  However,  the  focus  here  is  upon  linear 
equilibrium  sorption  to  promote  a  ready  and  intuitive  grasp  of  the  pertinent  con- 
cepts. For  this  case,  /i  may  be  thought  of  as  the  time  that  the  solute  is  sorbed.  Thus 
the  arrival  time  of  a  particle  undergoing  a  linear  equilibrium  sorption  process  is  the 
sum  of  the  nonreactive  travel  time  r  and  the  time  sorbed  [i.  The  transfer  function 
characterizing  this  transport  and  reaction  process  is  (see  Cvetkovic  et  al.  [1998]) 

j[t]x}  =  S[t-r[x}-fi[x}}  (4.3) 

Following  the  work  of  Cvetkovic  et  al.  [1998],  the  adopted  reaction  parameter 
model  is 

P  =  Pg  exp[w]  exp[/?  In  k]  (4.4) 

where  Pg  is  the  geometric  mean  of  P,  w  is  a  zero-mean,  exponentially-correlated 
normal  random  variable  with  variance  cr^  and  correlation  length  Xw  =  X\nk  that  is 
uncorrelated  to  the  log  conductivity  field,  and  )3  is  a  strength  of  correlation  parameter 
relating  P  to  the  underlying  log  conductivity  field.  The  equality  Xw  =  Ainfc  was 
specified  for  convenience. 


82 

Two  specific  injection  modes  are  considered:  uniform  injection  in-flux  and  uniform 
resident  injection.  For  the  uniform  injection  in-flux,  each  streamtube  receives  an  equal 
amount  of  solute  mass  (Equation  (3.5)),  and  thus  mass  weight  (Equation  (3.6)).  Uni- 
form resident  injection  distributes  the  mass  uniformly  in  space,  but  each  streamtube 
receives  a  different  mass  and  mass  weight  (Equation  (3.3)).  The  mass  breakthrough 
for  some  streamtube  i  is  given  by 

mi[t;x]=mi[0;0)j[t-x]  (4.5) 

An  ergodic  condition  is  assumed  such  that  aggregate  properties  may  be  thought  of 
as  ensemble  properties  in  the  statistical  sense.  The  aggregate  mass  breakthrough 
is  given  by  taking  the  expected  value  of  Equation  (4.5).  For  injection  in-flux,  the 
expected  solute  breakthrough  curve  is  given  by 

<Sf[t]>=  M  <7M>  (4-6) 

For  uniform  resident  injection,  the  expected  solute  breakthrough  curve  is  given  by 

<  Su[t]  >=  MU  <  7[t;  x]/v[0]  >  (4.7) 

where  U  is  the  arithmetic  mean,  or  spatially  averaged,  velocity  and  v[0]  is  the  velocity 
at  the  injection  point.  The  temporal  moment  of  this  breakthrough  curve  is  sought. 
The  n  temporal  moment  is  given  by 

tn  =  lim(-l)n|^<,SM>  (4.8) 

s->0  QSn 

where  S  is  the  Laplace  transform  of  S  and  s  the  corresponding  Laplace  space  variable. 
The  Laplace  transform  of  Equation  (4.3)  is  given  by 

j[s;  x]  =  exp  [-s  (r[x]  +  fi[x})}  (4.9) 
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Two  temporal  moments  of  interest  for  the  in-flux  breakthrough  curve  are  the  mean 

tlf  =<t>  +  <n>  (4.10) 

and  the  variance 

t2f  =  a2T  +  al  +  2aT^  (4.11) 
For  uniform  resident  injection,  the  moments  are 

tlf  =<  Ut/v[0]  >  +  <  Un/v[0]  >  (4.12) 

and  the  variance 

t2f  =<  U{t  +  n)2/v[0]  >  -  <  Ut/v[Q)  >2  -  <  Ufi/v[0]  >2  (4.13) 

From  these  equations,  it  is  seen  that  the  mean  and  variance  of  the  expected  reactive 
travel  time  distribution  is  comprised  of  the  first  two  moments  of  the  nonreactive  travel 
time  and  the  reaction  flow  path  and  the  cross  moment  of  the  travel  time  and  reaction 
flow  path. 

4.4    Simulation  of  Reactive  Solute  Transport 

The  results  of  the  simulations  described  in  Chapter  2  are  employed  to  simulate  the 
transport  of  a  reactive  solute.  As  mentioned  in  that  chapter,  information  was  recorded 
along  two  particle  trajectories:  the  first  was  a  natural  trajectory  corresponding  to  a 
hydrodynamic  streamline,  and  the  second  was  a  straight  line  parallel  to  the  mean 
flow  direction.  Data  were  recorded  at  equal  increments  in  time  after  injection,  and 
at  "control  planes"  perpendicular  to  the  mean  flow  at  equal  displacements  from  the 
injection  plane.  Among  these  data  were  the  velocity,  the  hydraulic  conductivity,  and 
an  observation  of  a  lognormal  random  field  with  the  same  statistical  and  correlation 
properties  as  the  log  hydraulic  conductivity  field. 
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4.4.1    Reaction  Parameter 

From  these  data,  trajectory-based  P  fields  were  constructed  using  Equation  (4.4). 
In  this  manner,  any  variety  of  correlation  relationships  may  be  analyzed  using  the 
same  set  of  simulation  data.  There  is  a  dearth  of  information  about  the  structure  of 
sorption  fields  in  natural  systems  (see,  e.g.  Tompson  [1993]  or  Jawitz  [1999]).  Thus, 
numerically  convenient  parameter  values  were  selected.  This  suite  of  experiments  ex- 
amined correlation  parameter  values  of  /3  =  {-1,  -1/2, 0, 1/2, 1}  in  log  conductivity 
fields  of  variances  ofnJt  =  {1/2, 1}.  The  extreme  correlation  parameters  {3  =  {—1,1} 
are  arguably  physically  implausible,  yet  they  are  included  to  demonstrate  the  effect 
of  correlation  without  the  "noise"  associated  with  an  uncorrelated  component.  The 
remaining  three  values,  f)  =  {—1/2, 0, 1/2}  might  be  considered  a  range  of  more  plau- 
sible values.  For  consistency,  the  variance  of  point  log  reaction  parameter  was  held 
constant  and  equal  to  that  of  the  log  conductivity  field  by  requiring  a2w  =  (1  — /32)ofnA.. 
The  expected  value  of  P  was  1.  This  is  a  slightly  different  approach  than  that  taken 
with  the  log  hydraulic  conductivity,  where  the  geometric  mean  was  fixed  at  Kg  =  1. 
In  the  case  of  the  two-dimensional  isotropic  log  conductivity  field,  the  effective  con- 
ductivity is  equal  to  the  geometric  mean.  The  expected  value  of  the  reaction  value 
was  fixed  <P>=  1,  because  in  some  practical  applications  P  is  assumed  to  be  lin- 
early related  to  some  property  of  interest,  and  the  "amount"  of  that  property  is  to 
be  estimated  based  upon  the  temporal  moments  of  linearly  sorbing  solutes,  or  tracers 
(see  e.g.,  Jin  et  al.  [1995],  Annable  et  al.  [1998],  and  Jawitz  [1999]). 

The  point  average  of  reaction  parameter  P  as  observed  along  a  Lagrangian  trajec- 
tory is  a  function  of  injection  mode,  correlation  to  the  log  conductivity  and  the  vari- 
ance of  the  log  conductivity.  As  with  the  log  conductivity,  the  equal-flux  streamtubes 
exhibit  an  apparently  stationary  mean  reaction  parameter  in  displacement  along  the 
mean  flow  direction,  whereas  the  equal-area  streamtubes  exhibit  distinctly  nonsta- 
tionary  behavior  that  is  analogous  to  that  of  the  log  conductivity  (see  Figures  4.1  and 
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4.2).  These  observations,  of  course,  are  not  surprising,  since  the  reaction  parameter  is 
a  function  of  the  log  conductivity.  Positive  correlation  between  the  log  conductivity 
and  the  reaction  parameter  (/3  >  0)  results  in  a  reaction  parameter  trajectory  mean 
that  is  larger  than  the  uniform  spatial  mean.  Negative  correlation  (/3  <  0)  results  in 
a  reaction  parameter  trajectory-mean  that  is  smaller  than  the  uniform  spatial  mean. 
No  correlation  (/?  =  0)  results  in  equal  trajectory-based  and  uniform  spatial  mean 
reaction  parameter  estimates.  A  decrease  in  the  system  variability  (e.g.,  the  variance 
of  log  conductivity  ofnk)  results  in  a  smaller  difference  between  the  volume  average 
and  the  trajectory  averages. 

4.4.2    Reaction  Flow  Path 

A  simple  trapezoidal  rule  was  used  to  calculate  the  reaction  flow  path  fi  along  a 
solute  particle  trajectory.  For  a  collection  of  n  trajectories,  there  were  n  \i  "obser- 
vations" at  each  control  plane.  For  an  intuitive  grasp  of  what  the  statistics  of  the 
reaction  path  "means" ,  recall  that  for  linear  equilibrium  sorption,  the  reaction  path 
corresponds  to  the  time  that  the  solute  is  sorbed.  The  sum  of  the  the  nonreactive 
travel  time  and  the  sorbed  time  equals  the  total  reactive  travel  time  of  the  solute 
from  the  point  of  injection  to  the  point  of  observation.  See  Cvetkovic  et  al.  [1998]  for 
further  application  of  the  statistics  of  the  reaction  flow  path. 

The  mean  and  the  variance  of  the  reaction  flow  path  are  presented  for  Lagrangian 
trajectories  that  carry  either  an  equal  flow  weight  or  an  equal  area  weight  at  a  ref- 
erence plane.  Recall  that  the  equal  flow  collection  can  be  thought  to  correspond 
to  a  uniform  injection  in-flux  and  that  the  equal  area  collection  can  be  thought  to 
correspond  to  a  uniform  injection  into  the  resident  fluid.  The  coupled  effect  of  injec- 
tion mode  and  reaction  parameter  correlation  to  the  log  conductivity  are  illustrated 
by  examining  reactive  flow  fields  with  five  different  correlation  characteristics.  The 
first  case  is  that  in  which  the  reaction  parameter  is  completely  determined  by  the 
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Figure  4.1:  Mean  reaction  parameter  values  for  different  values  of  correlation  pa- 
rameter P  and  different  trajectory  collections  for  the  afnk  =  1  set  of  simulations. 
Top  figure  is  /?  =  {-1,+1}.  Middle  figure  is  f3  =  {-1/2, +1/2}.  Bottom  figure  is 
/3  —  0  (no  correlation).  Trajectory  collections  are  area-weighted  Lagrangian  (AL), 
flow- weighted  Lagrangian  (AF)  and  area- weighted  Eulerian  (AE). 
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Figure  4.2:  Mean  reaction  parameter  values  for  different  values  of  correlation  pa- 
rameter (3  and  different  trajectory  collections  for  the  ofnfc  =  1/2  set  of  simulations. 
Top  figure  is  P  =  {-1,+1}.  Middle  figure  is  0  =  {-1/2, +1/2}.  Bottom  figure  is 
P  —  0  (no  correlation).  Trajectory  collections  are  area-weighted  Lagrangian  (AL), 
flow- weighted  Lagrangian  (AF)  and  area- weighted  Eulerian  (AE). 
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conductivity  field  (ft  =  {-1,1}).  Considered  are  the  two  sub-cases,  namely  that  of 
"perfect"  negative  correlation  (/5  =  -1)  and  "perfect"  positive  correlation  (/5  —  -1). 

For  /3  =  —  1,  the  equal  area  collection  exhibits  nonlinear  behavior  analogous  to 
that  of  travel  time  (see  Figure  4.3  and  4.4).  However,  this  nonlinearity  is  enhanced  by 
the  correlation  of  high  reaction  parameter  values  in  the  low  conductivity  areas.  For 
the  equal  flow  collection,  however,  the  nonlinearity  is  mitigated,  and  \i  appears  to  be 
a  linear  function  of  displacement.  For  /3  =  1,  the  both  trajectory  collections  exhibit 
the  same  linear  behavior.  Thus  positive  correlation  appears  to  mitigate  the  effect  of 
injection  mode  on  the  mean  behavior.  The  reaction  path  variance  is  strongly  affected 
by  the  sign  of  the  correlation  (see  Figures  4.5  and  4.6).  For  negative  correlation,  there 
is  a  distinct  difference  between  the  injection  mode.  However,  positive  correlation 
again  mitigates  the  difference  in  modes.  Moreover,  the  positive  correlation  greatly 
reduces  the  reaction  path  variance.  For  linear  equilibrium  sorption,  this  implies  a 
reduction  of  the  reactive  travel  time  variance. 

The  next  case  considered  is  that  of  no  correlation  between  the  reaction  param- 
eter and  the  log  conductivity  (/?  =  0).  The  reaction  path  mean  exhibits  the  same 
behavior  as  does  nonreactive  travel  time.  The  nonlinearity  in  the  mean  for  the  equal 
area  streamtubes  is  attributable  to  the  nonstationarity  in  the  associated  velocity  field. 
However,  the  combined  effects  of  heterogeneous  velocity  and  sorption  manifest  them- 
selves in  the  reaction  flow  path  variance.  The  reaction  flow  path  variance  is  higher 
than  that  of  the  nonreactive  travel  time,  even  though  the  scaled  means  are  the  same. 
The  /3  —  0  reaction  flow  path  variance  "lies  in  between"  that  of  the  positive  and 
negative  correlation.  For  the  lower  "total  system  variability"  case  of  afnk  =  1/2,  the 
magnitude  of  the  effect  of  correlation  appears  nearly  symmetric  (see  Figure  4.6).  For 
cinfc  —  15  however,  there  the  effect  of  negative  correlation  seems  to  have  a  greater 
impact  upon  both  total  variability  and  the  effect  of  injection  mode  (see  Figure  4.5). 
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Figure  4.3:  Mean  reaction  path  observed  values  and  estimators  for  different  values 
of  correlation  parameter  ft  and  different  trajectory  collections  for  the  crfnk  =  1  set  of 
simulations.  Top  figure  is  ft  =  {-1,  +1}.  Middle  figure  is  /3  =  {-1/2,  +1/2}.  Bottom 
figure  is  ft  =  0  (no  correlation).  Trajectory  collections  are  area- weighted  Lagrangian 
(AL)  and  flow- weighted  Lagrangian  (AF). 
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Figure  4.4:  Mean  reaction  path  observed  values  and  estimators  for  different  values  of 
correlation  parameter  and  different  trajectory  collections  for  the  afnk  =  1/2  set  of 
simulations.  Top  figure  is  P  =  {-1,+1}.  Middle  figure  is  /?  =  {-1/2,  +1/2}.  Bottom 
figure  is  /3  =  0  (no  correlation).  Trajectory  collections  are  area- weighted  Lagrangian 
(AL)  and  flow- weighted  Lagrangian  (AF). 
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Figure  4.5:  Reaction  path  variance  observed  values  and  estimators  for  different  values 
of  correlation  parameter  @  and  different  trajectory  collections  for  the  afnk  =  1  set  of 
simulations.  Top  figure  is  /3  =  {— 1,+1}.  Bottom  figure  is  /3  =  {—1/2, +1/2}.  Tra- 
jectory collections  are  area-weighted  Lagrangian  (AL)  and  flow-weighted  Lagrangian 
(AF). 
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Figure  4.6:  Reaction  path  variance  values  and  estimators  for  different  values  of  cor- 
relation parameter  /3  and  different  trajectory  collections  for  the  afnk  =  1/2  set  of 
simulations.  Top  figure  is  /3  =  {-1,+1}.  Bottom  figure  is  /3  =  {-1/2,-1-1/2}.  Tra- 
jectory collections  are  area-weighted  Lagrangian  (AL)  and  flow-weighted  Lagrangian 
(AF). 
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4.4.3    Reaction  Flow  Path-Travel  Time  Cross  Moment 

The  reaction  flow  path-nonreactive  travel  time  cross  moment  a^T  contributes  to 
the  variance  of  the  expected  solute  flux  (see  Equations  (4.11)  and  (4.13)).  For  linear 
equilibrium  sorption,  this  reflects  the  correlation  between  the  time  sorbed  and  the 
time  traveling.  The  correlation  of  the  reaction  flow  path  and  the  travel  time  is  little 
affected  by  injection  mode  (see  Figures  4.7  and  4.8).  It  is,  however,  sensitive  to  the 
correlation  parameter  @.  As  with  the  reaction  flow  path  statistics,  the  deviation  of 

from  the  /?  =  0  is  greater  in  the  case  of  positive  correlation  for  the  same  value 
of  /?.  The  "total  system  variability"  seems  to  diminish  the  effect  of  correlation  upon 
aMT.  This  is  most  readily  apparent  from  a  comparison  of  the  two  ft  =  {—1,1}  cases 
for  the  different  values  of  afnk  (compare  the  top  figures  in  Figures  4.7  and  4.8).  The 
correlations  for  the  correlated  case  "shift"  towards  that  of  the  /3  —  0  in  the  higher  log 
conductivity  variance  field,  and  the  f3  =  0  case  is  slightly  lower  for  the  greater  afnk.  In 
all  cases  except  /3  =  l,<jfak  =  1,  the  correlations  appear  to  plateau  at  nonzero  values. 
This  is  not  surprising,  since  both  r  and  n  are  integrals  of  the  inverse  Lagrangian 
velocity. 

4.5    Theory  of  Reactive  Solute  Transport 

The  concepts  of  the  previous  chapters  are  applied  to  extend  the  work  of  researchers 
including  Cvetkovic  et  al.  [1998]  and  Demmy  et  al.  [1999].  These  results  are  used 
to  interpret  the  results  of  the  numerical  experiments.  Expressions  are  sought  for 
the  mean  and  variance  of  the  reaction  flow  path  for  the  two  different  collections 
of  streamtubes.  Cvetkovic  et  al.  [1998]  have  developed  expressions  for  equal  area 
streamtubes  in  unbounded  multidimensional  flow  fields.  Their  mean  \i  expression 
neglects  the  correlation  between  the  reaction  parameter  and  the  log  conductivity, 
and  is  a  linear  function  of  displacement.  Their  variance  expression  is  based  upon 
a  series  expansion  of  the  inverse  Eulerian  velocity  perturbation.  This  work  extends 
this  work  with  consideration  the  correlation  between  the  reaction  parameter  and  log 
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Figure  4.7:  Correlation  between  the  reaction  flow  path  and  travel  time  for  different 
values  of  correlation  parameter  /?  and  different  trajectory  collections  for  the  ofnJfc  =  1 
set  of  simulations.  Top  figure  is  j3  =  {-1,+1}.  Bottom  figure  is  (3  =  {-1/2, +1/2}. 
Trajectory  collections  are  area-weighted  Lagrangian  (AL)  and  flow-weighted  La- 
grangian  (AF). 
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Figure  4.8:  Correlation  between  the  reaction  flow  path  and  travel  time  for  differ- 
ent values  of  correlation  parameter  /3  and  different  trajectory  collections  for  the 
a\nk  —  1/2  set  of  simulations.  Top  figure  is  /3  =  {— 1,+1}.  Bottom  figure  is 
ft  —  {—1/2, +1/2}.  Trajectory  collections  are  area-weighted  Lagrangian  (AL)  and 
flow- weighted  Lagrangian  (AF). 
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conductivity  in  the  calculation  of  the  mean,  working  from  Lagrangian  properties  that 
are  stationary  in  displacements  along  the  mean  flow  direction,  and  examination  of 
the  effects  of  injection  mode. 

4.5.1    Mean  Reaction  Flow  Path 

Substitution  of  Equation  (4.4)  into  Equation  (4.2)  and  taking  the  expected  value 
yields 


< 


The  stationary  Lagrangian  serves  as  a  reference  field.  The  log  conductivity  field  is 
decomposed  into  its  mean  and  mean  removed  perturbation  /  and  the  velocity  into 
its  mean  and  mean  removed  perturbation  thus 


< 


Define  Pg  =  P*  exp[/xinfc],  and  expand  the  nonlinear  terms  of  the  previous  equation 
into  series,  take  the  expected  value  retaining  second  order  terms  thus 

<M>=  P*exp^/2](1  -  fa,  +  fo}/2  +  al)x  (4.16) 

Recall  that  aj  =  afnk.  The  /  subscript  is  used  for  clarity.  This  expression  captures 
the  linearity  of  the  observed  mean  \x  (see  Figure  4.3),  but  shows  a  heavy  dependence 
upon  the  sign  of  the  correlation  parameter  f3  that  is  not  apparent  in  the  data.  The 
mean  reaction  flow  path  for  the  equal  flow  streamtubes  appears  to  be  dependent  only 
upon  the  mean  value  of  the  reaction  parameter,  e.g. 

<»>=P>^Kll+^>l\+<Z)X  (4,7) 
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For  the  hypothesized  Lagrangian  distribution  V  =  U(l  +  baj)  and  a\  =  baj,  thus  the 
Equation  reduces  to 

<»>=  PgexV[a2j2  +  P2a2/2]x/U  (4.18) 
The  equal  area  streamtube  mean  [i  is 

<W„[„]>=p;eXp[4/2,,f(2^)^  (4,9) 
Manipulations  similar  to  that  which  yielded  Equation  (4.16)  yield 

<  Uv/v[0]  >=  ^^Hf  f  (1  +  2a„2  +  P2o2/2  +  Cv[x]  -  2(3Cvf[x})dx  (4.20) 

1/(1+00"  j)  JQ 

Substituting  the  simplified  relationship  Cvf[x]  =  a2  exp[-a;/(2A)]/2  yields 

<w»[°i>40^<1+2^+^2)*  (4.21) 

+  2/3AaJ(exp[-x/(2A)]  -  1)  +  \o){\  -  exp[-ox/A]) 
Notice  that  the  first  collection  of  terms,  or  the  linear  portion,  of  the  preceding  equa- 
tion is  the  same  as  that  in  Equation  (4.16).  Similar  arguments  are  employed  to  yield 
the  equal  area  streamtube  mean  //  estimator 

Pgexp[^/2  +  /?V;/2] 
<  U,/v[0]  >=  f—x  (4  22) 

+  2^AaJ(exp[-x/(2A)]  -  1)  +  Xa){\  -  exp[-6x/A]) 

4.5.2    Reaction  Flow  Path  Variance 

Adapting  an  expression  for  the  reaction  flow  path  variance  given  by  Cvetkovic 
et  al.  [1998]  (their  Equation  (B4))  to  the  simplified  velocity  and  velocity-log  conduc- 
tivity functions  central  to  this  work,  write 

al\x]  =  2P2  exp[a2w}/U2  /  (x  -  s)(Cw[s]  +  p2Cf[s]  -  4pCfv[s]  +  Cv[s})ds  (4.23) 

Jo 
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For  the  special  case  of  a\  = 


(4.24) 


+  (2-4P)x/\f} 


A  simple  and  intuitive  reaction  flow  path  variance  estimator  for  the  equal-area 
streamtube  collection  is  untenable  in  this  framework.  This  is  due  in  part  to  the 
nonlinear  reaction  flow  path  mean. 


Figures  4.3  and  4.4  demonstrate  a  close  match  between  Equations  (4.18)  and 
(4.22)  and  simulation  data.  The  parameters  placed  into  the  estimators  are  the  target 
values  selected  for  the  simulation,  and  thus  represent  an  a  priori  estimate  of  the 
mean  reaction  flow  path.  Equation  (4.22)  captures  not  only  the  nonlinearity  of  the 
nonreactive  travel  time,  but  the  nonlinearities  induced  by  reaction  parameter  and  log 
conductivity  correlation.  While  the  "absolute  difference"  between  Equations  (4.18) 
and  (4.22)  appear  small,  consider  that  the  scale  of  the  simulations  is  fairly  large.  The 
solute  has  traveled  over  forty  correlation  scales  of  the  log  conductivity.  In  settings 
where  the  the  scale  of  conductivity  heterogeneity  are  large  compared  to  the  overall 
scale,  the  relative  impact  of  injection  mode  and  reaction  parameter  correlation  to  the 
log  conductivity  will  have  a  relatively  greater  effect. 

The  equal-flux  streamtube  reaction  flow  path  variance  estimator  (Equation  (4.24)) 
shows  relatively  poor  performance  in  the  more  heterogeneous  system  (Figure  4.5), 
compared  to  that  in  the  less  heterogeneous  system  (Figure  4.6).  Moreover,  this  ex- 
pression does  not  work  well  at  all  for  values  of  positive  correlation,  yielding  physically 
implausible  reaction  flow  path  variances  (cr^  <  0)  for  correlation  parameter  values 
greater  that  1/2.  The  lack  of  "robustness"  is  perhaps  attributable  to  an  inadequate 


4.6  Discussion 
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representation  of  the  myriad  interacting  processes  by  the  first-order  perturbation  ex- 
pansions. 

Injection  of  solute  in  the  influent  fluid  mitigates  the  effect  of  correlation  on  the 
mean  and  the  nonlinearity  of  mean  nonreactive  travel  time.  For  tracer  experiments 
that  determine  the  volume  or  area  of  some  property  based  upon  mean  arrival  times 
of  reactive  and  nonreactive  tracers,  correlation  will  have  little  effect  upon  the  first 
temporal  moment  of  the  solute  breakthrough  curve  if  and  only  if  the  solute  is  uni- 
formly injected  in  flux.  However,  injections  that  preferentially  favor  low-flux  areas 
may  exhibit  strong  nonlinearities,  and  be  particularly  sensitive  to  negative  correla- 
tions between  the  local  conductivity  and  reaction  parameter. 

The  contribution  of  the  fj,  variance  to  the  second  moment  of  the  reactive  solute 
travel  time  is  extremely  sensitive  to  correlation  between  the  reaction  parameter  and 
the  conductivity.  Experiments  that  rely  upon  higher  order  moments  of  tracer  ex- 
periments are,  therefore,  questionable  if  there  is  no  p-ln  k  correlation  information. 
Whereas  injection  mode  dominates  the  correlation  behavior  in  the  mean,  correlation 
dominates  injection  mode  in  the  reactive  travel  time  variance. 

4.7    Spatial  Moments 

Symmetry  would  dictate  proceeding  with  a  corresponding  spatial  moment  analysis 
for  transport  of  solute  subject  to  a  linear  equilibrium  sorption  in  a  heterogeneous 
reaction  parameter  field.  However,  this  process  presents  a  particular  complication 
that  thwarts  the  development  of  relatively  simple  expressions  following  the  same 
methods  those  previously  employed.  In  the  case  of  the  travel  time  statistics,  the 
"control  structures"  are  two  control  planes.  A  fluid  and  solute  parcel  travels  between 
these  two  control  planes,  sweeping  out  a  fixed  and  constant  volume,  without  regard  to 
correlation  properties,  reaction  parameters,  or  processes.  The  time  required  to  sweep 
this  volume,  of  course,  is  highly  dependent  upon  these  factors.  For  a  "clock  time" 
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control  structure,  a  solute  parcel  will  trace  different  volumes  for  different  processes 
and  correlation  structures. 

The  reaction  parameter  model  is  comprised  of  two  random  processes.  The  first 
is  that  of  the  log  conductivity  field  and  the  second  is  a  random  component  that 
has  the  same  distributional  properties  as,  but  is  otherwise  uncorrected  to,  the  log 
conductivity  field.  In  the  case  of  the  temporal  moments,  the  weighted  covariance 
functions  of  the  log  conductivity  and  this  random  component  were  employed  directly. 
These  functions  can  be  parameterized,  of  course,  in  terms  of  the  nonreactive  travel 
time.  However,  these  travel  time  parameterized  covariance  functions  time  are  not 
immediately  applicable  to  sorption  analysis  in  the  way  that  covariance  parameterized 
in  clock  time  would  be,  because  the  displacement  for  a  given  clock  time  is  highly 
dependent  upon  the  correlation  of  the  reaction  parameter  to  the  log  conductivity 
field.  This  is  not  to  say  that  sorption-dependent  spatial  moment  analysis  is  not  viable 
in  this  framework.   In  fact,  Cvetkovic  et  al.  [1998]  analyzed  the  spatial  moments 
of  a  solute  subject  to  nonlinear  equilibrium  sorption  in  a  heterogeneous  sorption 
field  in  this  framework.   However,  the  sorption  parameter  and  the  underlying  log 
conductivity  field  were  uncorrelated.  Further  clock  time  oriented  are  developments 
left  for  future  work.  Expressions  involving  the  travel  time  dependent  reaction  flow 
path  are  developed  and  applications  to  and  implications  for  first  order  decay  are 
discussed. 

Continuing  the  theme  of  this  dissertation,  expressions  based  upon  the  trajectory- 
based  covariance  functions  are  sought.  The  time-dependent  reaction  flow  path  for 
the  selected  reaction  parameter  model  is  related  to  the  log  conductivity  covariance 
function,  as  shall  be  demonstrated  in  the  following. 

As  in  the  previous  chapter,  transport  along  Lagrangian  trajectories  is  considered. 
Assumed  are  second-order  stationary  velocities  along  area-weighted  trajectories  in 
time  and  second-order  stationary  velocities  along  flow-weighted  trajectories  in  space. 
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The  log  hydraulic  conductivity  is  assumed  to  be  similarly  stationary.  Moreover,  the 
point  random  process  Pg  exp[w]  is  assumed  likewise  stationary,  since  it  is  uncorrelated 
to  the  random  log  conductivity.  From  these  assumptions,  the  expected  value  of  the 
time-dependent  reaction  flow  path  is 

<  h[t]  >  =  ^  <  Pg  exp[w[t}]  exp[p  In  k[t}}  >  dt 
=  Pgexp[a2j2  +  /32alk/2}r 

The  variance  is  given  by 

var[/x[r]]  =<  f  Pg  exp  [w[t]]  exp  [0  In  k[t}}  dt 

Jo  (4-26) 

/  Pg  exp  [w[t\]  exp  [P  In  k[t]}  dt'-  <  h[t]  >2 
J  o 

Consider  the  integrals  in  the  previous  equation.  The  exact  expression  is  denoted  by 
I  and  is  approximated  by  the  truncated  series  expansions  of  the  exponential  terms 
thus 

J~pg  f  [T  <  {I +  w[t})  {I +  w[t'])  (1+P\n  k[t})  (1+ 13  In  k[t'])>dtdt'  (4.27) 
Jo  Jo 

Note  that  a  possibly  more  elegant  approach  would  be  to  explore  the  Lagrangian 
features  of  the  conductivity  field,  rather  than  the  log  conductivity  field,  and  work 
from  a  stationary  conductivity  field  covariance  function  and  a  similar  uncorrelated 
field.  This  would  obviate  the  need  for  the  expansions.  This  suggested  inquiry  is  left 
to  future  works.  Expanding  the  above  expression,  taking  the  expected  value,  and 
dropping  products  of  covariances  yields 

I  =  P2g  f  f  (1  +  Cw[t,t'}  +  (32Clnk[t,t'})dtdt'  (4.28) 
7o  Jo 

A  few  observations  allow  the  adoption  of  an  approximate  expression  for  the  time- 
dependent  Lagrangian  covariance  function  C\nk[t,  t'].  Relatively  simple  expressions 
that  capture  the  behavior  of  the  dominant  processes  are  sought  here.  These  expres- 
sions might  facilitate  future  works  that  deemphasize  simplicity  for  greater  rigor  and 
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perhaps  a  better  fit  to  the  data.  The  Eulerian  log  conductivity  covariance  is  the 
foundation  upon  which  the  rest  of  this  work  has  been  built.  The  model  field  is  ex- 
ponentially correlated  in  space,  with  correlation  length  Ainfc.  The  Eulerian  velocity 
covariance  exhibits  a  much  longer  correlation  length  due  to  the  dissipation  of  the 
head  across  the  field  (see  previous  chapters).  It  is  reasonable,  therefore,  to  expect 
that  the  Lagrangian  log  conductivity  covariance  correlation  length  will  lie  somewhere 
in  between  that  of  the  Eulerian  log  conductivity  and  the  Lagrangian  velocity. 

At  zero  lag,  the  area-weighted  Lagrangian-trajectory  log  conductivity  covariance 
is  equal  to  the  variance  of  the  Eulerian  field.  This  value  is  afnfc,  by  definition.  Thus, 
an  exponential  approximation  of  the  Lagrangian  log  conductivity  covariance  function 
is 

ClDkiat[t]  =  alk  exv[-b'Ut/\lnk}  (4.29) 

where  b'  is  a  function  similar  to  b  that  relates  the  correlation  length  to  a  correla- 
tion time.  The  approximation  b'  —  1  is  made,  which  is  equivalent  to  assuming  that 
the  Lagrangian  and  Eulerian  correlation  lengths  are  equal  is  space  and  substituting 
x  =  Vht  =  Ut  into  the  space-dependent  Eulerian  log  conductivity  covariance  func- 
tion. Comparison  of  Equation  (4.29)  with  b'  =  1  to  covariance  functions  estimated 
from  log  conductivity  observations  taken  along  area-weighted  Lagrangian-trajectories 
generated  in  the  numerical  experiments  described  in  Chapter  2  shows  good  qualita- 
tive agreement  (see  Figure  4.9).  The  data  exhibit  an  apparent  dependence  of  the 
Lagrangian  correlation  length  upon  the  variance  of  the  log  conductivity  field  that  is 
not  captured  by  Equation (4. 29).  One  might  expect  similar  behavior  in  the  Eulerian 
field  a  priori,  since  the  appropriate  substitution  into  the  Eulerian  covariance  function 
would  be  x  =  Uht  =  Ut/(1  +  bafnk).  The  higher  variance  of  log  conductivity  results 
in  a  lower  harmonic  mean  average  velocity,  implying  it  takes  longer  to  to  travel  a 
unit  length  in  more  heterogeneous  fields  and  a  larger  Lagrangian  correlation  time.  In 
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Figure  4.9:  Assumed  time-dependent  log  conductivity  covariance  function  along  La- 
grangian trajectories  (Equation  (4.29))  compared  to  simulation  data  for  variances 
of  log  conductivity  afnk  =  {1/2,1}.  The  Lagrangian  velocity  covariance  function  is 
Equation  (2.61). 

fact,  this  is  the  behavior  observed  in  the  simulations  (see  Figure  4.10)  The  data  are 
compelling  enough  to  suggest  further  work.  However,  we  shall  accept  Equation  (4.29) 
with  b'  =  1  as  a  suitable  approximation. 

Substituting  Equation  (4.29)  with  t  =  \t  -  t'\  into  Equation  (4.28),  integrating 
and  inserting  the  result  into  Equation  (4.26)  yields 


var[/i[r]]  = 


U2 


iUT/Kk-(l-exp[-Ur/Xlnk}) 
-P?exp[al  +  p2olk]T2 


(4.30) 


Neither  this  equation  nor  Equation  (4.25)  show  a  dependence  upon  the  sign  of  (3. 

For  the  mean  reaction  flow  path,  the  simulations  show  a  weak  dependence  upon 
correlation  parameter  (see  Figure  4.11).  The  linear  propagation  of  the  mean  reaction 
flow  path  in  time  is  due  to  the  stationarity  of  the  local  reaction  parameter  along 
the  area-weighted  Lagrangian  trajectories.  Again,  it  should  be  noted  that  "time" 
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Figure  4.10:  Hypothesized  time-dependent  Eulerian  log  conductivity  covariance  func- 
tions for  ofnJfc  =  {1/2, 1}  compared  to  simulation  data  and  Equation  (4.29). 


Figure  4.11:  Mean  reaction  flow  path  as  a  function  of  nonreactive  travel  time  com- 
pared to  estimator  Equation  (4.25)  for  correlation  parameter  values  /5  =  {  —  1,0, 1}. 
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refers  to  nonreactive  travel  time.  Reaction  flow  path  variance,  however,  is  strongly 
affected  by  correlation,  and  Equation  (4.30)  is  completely  inappropriate  (see  Fig- 
ure 4.12).  Moreover,  this  equation  can  give  physically  implausible  values  at  large 
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Figure  4.12:  Reaction  flow  path  variance  a  function  of  nonreactive  travel  time  for 
correlation  parameter  values  /3  =  {—1,0, 1}. 

values  of  time.  This  is  a  result  of  the  series  approximations  required  for  the  "raw 
second  moment"  term  in  the  variance,  and  the  lack  of  such  approximations  in  the 
"first  moment  squared"  term. 

4.7.1    Injection  Mode 

The  coupled  effects  of  injection  mode  and  correlation  on  the  mean  travel  time 
dependent  reaction  flow  path  are  examined.  Using  the  relationships  developed  in 
Chapter  3,  the  hypothesis  that  the  statistics  reaction  flow  path  corresponding  to  the 
equal-flow  streamtubes  are  given  by  flow-weighting  the  equal-area  properties  is  made. 
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Thus,  the  mean  reaction  flow  path  for  a  flux  injection  is  given  by 

<ufrtT)>=<uf[pm>  (4.3i) 

In  taking  the  expected  value,  the  cross-correlation  of  the  initial  velocity  and  the  that 
portion  of  the  reaction  parameter  correlated  to  the  log  conductivity  field  must  be 
considered.  Using  now-familiar  expansions,  the  expected  value  of  the  flow- weighted 
reaction  flow  path  is 

<  ^W]  >=  ^^^^V  +  PCfu[t]/U)dt  (4.32) 
U  Jo 

where  Cfu  is  the  covariance  of  the  initial  velocity  perturbation  and  the  log  con- 
ductivity perturbation.  Once  again,  this  covariance  shall  be  approximated  with  an 
exponential  function  of  the  form  (see  Appendix  B) 

Ut 


(4.33) 


Substituting  this  term  into  Equation  (4.32)  and  integrating  yields 

<  f  ,H  >=  5^^!  (ur  +  ^.expI-J^])  (4.34) 


4.7.2  Discussion 

Similar  to  its  space  dependent  counterpart,  the  mean  time  dependent  reaction  flow 
path  is  sensitive  to  injection  mode.  However,  it  is  the  uniform  resident  injection  that 
appears  to  mitigate  correlation  effects  upon  the  mean  reaction  flow  path.  Correlation 
has  a  strong  effect  of  [i  variance,  and  it  is  quite  likely  that  these  effects  dominate  any 
injection  mode  effects,  as  was  the  case  for  the  space  dependent  \x. 
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Figure  4.13:  Mean  reaction  flow  path  as  a  function  of  nonreactive  travel  time 
compared  to  estimator  Equation  (4.34)  for  correlation  parameter  values  /3  = 
{-1/2,0,1/2}  (top  figure)  and  0  =  {-1,0,1}  (bottom  figure).  afDk  =  1. 


CHAPTER  5 
APPLICATION  OF  CONCEPTS 

5.1  Introduction 

An  illustration  of  a  practical  application  of  trajectory  based  flow  and  transport 
modeling  is  sought.  Useful  and  fairly  robust  relationships,  namely  predictors  of  spatial 
moments  of  a  continuously  injected  plume  subject  to  a  first  order  decay,  using  simpler 
streamtube  concepts  (e.g.  Jury  and  Roth  [1990]).  However,  the  observed  bias  in 
prediction  shall  be  explained  using  the  more  complete  trajectory  based  approach. 

In  light  of  the  robustness  of  the  simple  relationships  developed  here,  they  are  used 
to  estimate  relative  degradation  rates  for  different  solute  species  and  applied  to  field 
data  collected  at  a  field-scale  macrodispersion  and  natural  attenuation  experiment 
conducted  in  Columbus,  Mississippi  (see  Stauffer  et  al.  [1994]  and  Stauffer  et  al. 
[1997]). 

5.2  Theory 

Consider  a  subsurface  source  of  groundwater  contamination  that  releases  solute 
into  the  aqueous  phase  at  a  constant  concentration.  An  example  of  such  a  source  is 
a  poorly  mobile  non-aqueous  phase  liquid,  comprised  of  relatively  insoluble  organic 
compounds,  distributed  at  some  near-residual  saturation  through  out  a  portion  of 
some  aquifer.  The  "constant"  assumption  applies  when  the  rate  of  change  of  a  prop- 
erty of  interest,  such  as  the  aqueous  phase  concentration  at  the  source,  is  small  in 
comparison  to  the  temporal  scales  associated  with  transport  and  decay. 

The  actual  groundwater  flow  field  is  modeld  as  a  steady  irrotational  flow  field  that 
may  be  resolved  into  an  aggregation  of  streamtubes.  Neglecting  the  effect  of  local 
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dispersive  processes,  an  advection-decay  equation  with  streamtube  specific  parame- 
ters is  assumed  to  describe  the  fate  and  transport  of  the  solute  along  the  streamtube 
trajectory.  This  advection-decay  equation  is 

S-*  ™ 

subject  to  c[x,  t  <  0]  =  0  and  c[x  =  0,  t  =>  0]  =  Cq,  where  v  is  an  effective  advective 
velocity  [LT_1],  and  A:  is  a  first-order  decay  constant  [T-1].  The  assumption  of 
effective  streamtube  parameters  is  the  simplification  that  separates  this  method  from 
the  more  general  consideration  of  properties  varying  along  the  trajectory. 

Although  much  of  this  language  is  similar  to  that  used  in  earlier  chapters,  the 
fundamental  difference  is  that  the  streamtube  has  constant  effective  parameters.  This 
is  analgous  to  an  assuption  of  stationarity  of  properties,  but  is  different  in  significant 
ways,  as  shall  be  shown. 

For  advective  transport  in  a  single  streamtube,  the  resident  concentration  is  equal 
to  the  flux-averaged  concentration.  However,  the  difference  between  these  detection 
modes  manifests  itself  when  concentrations  are  averaged  over  several  streamtubes. 
For  the  remainder  of  this  discussion,  any  reference  to  a  concentration  will  implicitly 
imply  a  reference  to  a  resident  concentration.  The  solution  to  this  equation  for  the 
given  initial  and  boundary  conditions  is 

—  =  exp  [— kx/v]  h[vt  —  x]  (5.2) 

Co 

where  h[x]  is  the  Heaviside  step  function,  with  the  properties  h[x]  =  0  for  x  <  0  and 
h[x]  =  1  for  x  >  0. 

5.2.1    Spatial  Moments 

The  defintion  of  the  raw  spatial  moment  of  order  n  of  curve  c  is 

fj!n  =  J  oo°°xnc[a;](ix  (5.3) 
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A  normalized  moment  of  order  n  is  a  raw  moment  of  order  n  divided  by  the  zeroth 
moment,  and  is  denoted  by  dropping  the  prime  (i.e.,  /in).  The  zeroth  spatial  moment 
of  this  concentration  curve  is  given  by 

^[t]  =  l[l  -  exp[-fct]]  (5.4) 

The  zeroth  spatial  moment  of  a  contaminant  plume  related  to  the  total  mass  in 
solution.  The  first  normalized  spatial  moment,  or  center  of  mass,  is  given  by 

rtn  1  -  exp[-kt](kt  +  1) 
™  =  1*  J       1  -  exp[-!fct]  (  ] 

In  the  limit  t  — >  oo  n\  —  v/k. 

The  concentration  given  in  Equation  (5.2)  is  a  function  of  this  asymptotic  center 
of  mass  v/k.  Solving  Equation  (5.2)  with  h[-]  =  1  for  v/k  yields 

(5.6) 


k        In  [c/co] 

Thus,  an  observation  of  the  concentration  at  any  point  of  a  steady  streamtube  con- 
centration curve  is  an  observation  of  its  center  of  mass,  given  that  cq  and  x  are 
known. 

Assume  that  independent  samples  can  be  drawn  uniformly  from  the  collection  of 
streamtubes  that  originate  from  a  source  Cq  located  at  x  —  0.  From  the  linearity  of 
first  moments,  the  mean  center  of  mass  may  be  estimated  from  N  samples  by  the 
simple  estimator 

1    N  x- 

*"|>mT^  (5'?) 

5.2.2    Relative  Degradation 

Consider  two  solutes  a  and  b  subject  to  a  first-order  decay  process,  characterized 
by  rate  constants  A;a  and  kb,  respectively.  The  relative  degradation  rate  of  solute  a  to 
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b  is  defined  as 

K  =  ka/kb  (5.8) 

Let  these  solutes  be  released  into  a  flow  field  at  a  steady  rate,  and  allow  the  solutes 
to  reach  a  "steady  state"  in  which  the  concentration  profiles  for  each  do  not  vary  in 
time.  If  within  a  streamtube  i,  the  streamtube  effective  velocity  Vi  is  a  approximately 
equal  for  solutes  a  and  b,  the  relative  degradation  rate  of  a  to  b  in  streamtube  i  is 
given  by 

kri  =  (vi/kb)/(vi/ka) 

(5.9) 

That  is,  the  relative  degradation  rate  is  given  by  the  ratio  of  the  first  normalized 
spatial  moment  of  the  steady-state  concentration  curve  of  solute  b  to  that  of  a.  The 
average  relative  degradation  rate  is  estimated  as  the  arithmetic  mean  of  kri.  Notice 
that  this  is  not  the  ratio  of  the  mean  first  spatial  moments. 

5.3    Application:  Simulation 

The  estimator  was  tested  by  simulating  the  release  of  a  solute  into  a  heterogeneous 
two-dimensional  aquifer  using  the  United  States  Geological  Survey  flow  and  transport 
code  MOC3D  [Konikow  et  al,  1996].  The  length  units  used  herein  are  a  dimensionless 
quantity  resulting  from  the  normalization  of  the  fundamental  length  by  that  of  the 
correlation  length  of  the  exponentially  correlated  log  hydraulic  conductivity  field.  In 
the  interest  of  clarity,  however,  this  dimensionless  quantity  will  be  referred  to  as  a 
correlation  length  A. 

Let  the  reference  Cartesian  coordinate  system  be  oriented  such  that  the  ^-direction 
is  parallel  with  the  mean  direction  of  flow,  and  the  y-direction  is  orthogonal  to  that 
and  lies  within  the  plane  formed  by  the  aquifer.  Let  the  origin  lie  at  the  lower  left 
corner  of  the  aquifer.  The  rectangular  aquifer  extended  15  A  x  and  25  A  in  y,  and  was 
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subdivided  into  10  nodes  per  A  in  both  direction.  Constant  heads  maintained  at  the 
x  =  OA  and  x  =  15A  boundaries  imposed  a  steady  mean  gradient  of  J  =  0.01.  No-flow 
conditions  were  maintained  at  the  y  =  OA  and  y  =  25A  boundaries.  To  simulate  the 
release  of  a  contaminant  at  constant  concentration,  the  influent  water  concentration 
was  assumed  to  be  a  constant  concentration  c0  =  1.  This  "source"  was  centered  on 
the  x  =  0  boundary,  and  was  10A  wide. 

A  turning-bands  algorithm  generated  a  synthetic  heterogeneous  log  conductivity 
field  used  in  the  simulations  [Tompson  et  ai,  1989].  The  output  of  the  turning-bands 
generator  was  a  realization  of  an  exponentially-correlated  standard  normal  process 
(/x  =  0,cr  =  1).  This  field  was  subsequently  converted  to  fields  with  variances  of  0.1, 
0.5,  and  1.0  via  the  transformation 

/  =  VlnkZ  (5.10) 

where  a\a  *  is  the  target  standard  deviation  of  the  log  conductivity  field  /  and  z  is  the 
standard  normal  process.  The  aquifer  porosity  was  assumed  to  be  equal  to  the  mobile 
water  content  and  equal  to  a  constant  value  n  =  0.2.  The  effective  conductivity 
of  a  two-dimensional  exponentially  correlated  log  conductivity  field  is,  to  a  good 
approximation,  the  geometric  mean  conductivity  Kg.  For  \x  —  0,  Kg  =  1,  in  units  of 
(dimensionless)  correlation  lengths  per  characteristic  time.  The  independence  of  the 
effective  conductivity  from  its  variance  convieniently  allows  the  variance  to  be  changed 
without  substantially  changing  the  bulk  flow  properties.  Thus,  the  anticipated  Darcy 
flux  for  the  aquifer,  regardless  of  o\nk  was  =  0.01A  per  characteristic  time,  and  the 
characteristic  filtration  velocity  was  vi,  =  0.05A  per  characteristic  time.  The  subscript 
b  refers  to  a  bulk  property,  to  be  distinguished  from  streamtube  properties  discussed 
earlier. 

The  solute  was  assumed  to  undergo  a  constant  decay  process  characterized  by  a 
first-order  rate  coefficient  k  —  0.05  in  units  of  an  inverse  characteristic  time  defined 
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by  Vf/X.  This  value  yields  what  might  be  termed  a  "bulk"  first  spatial  moment 
Vb/k  =  1.  For  a  homogeneous  conductivity  field,  this  is  numerically  equivalent  to  the 
first  longitunial  spatial  moment  of  the  steady-state  plume. 

The  steady-state  concentration  field  was  "sampled"  by  taking  the  reported  nu- 
merical concentration  values  as  the  local  resident  concentration.  The  plume  center- 
of-mass  was  estimated  by  the  estimator  presented  in  Equation  (5.7).  The  "actual" 
first  longitudinal  spatial  moment  was  calculated  using  a  trapezoidal  rule  integration. 
A  fate  and  transport  simulations  were  carried  out  in  6  sets  of  hydraulic  conductivity 
fields  corresponding  to  6  standard  normal  process  realizations. 

5.3.1  Results 

For  the  case  of  a  homogeneous  conductivity  field,  any  one  non-zero  concentration 
measurement  used  with  the  estimator  given  in  Equation  (5.7)  returned  a  value  for  the 
first  longitudinal  spatial  moment  that  for  all  practical  purposes  was  equal  to  Vb/k,  and 
five  percent  less  than  that  predicted  by  the  numerical  integration.  As  the  variance  of 
log  conductivity  increased,  however,  both  the  numerical  integration  and  the  estimator 
tended  to  systematically  predict  values  greater  than  Vb/k.  Additionally,  the  estimator 
systematically  underpredicted  the  value  given  by  the  numerical  integration,  and  by 
an  amount  which  generally  increased  with  variance  of  the  log  conductivity. 

5.3.2  Discussion 

This  systematic  increase  in  the  the  first  longitudinal  spatial  moment  of  the  plume 
with  increasing  log  hydraulic  condutivity  variance,  despite  a  fairly  constant  vb  is 
explained  by  considering  the  nature  of  the  underlying  Lagrianian  velocity  field.  The 
steady-state  concentration  plume  is,  in  essence,  a  map  of  travel  times  from  the  source. 
As  previously  stated,  Cvetkovic  et  al.  [1996]  demonstrated  that  the  statisics  of  such 
a  field  are  nonstationary.  For  small  displacements  from  the  source,  the  mean  travel 
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Figure  5.1:  Resident  concentration  profile  averaged  along  a  transect  orthogonal  to 
the  mean  flow  direction  for  plumes  for  cr,2nfc  =  {0.0, 0.1, 0.5, 1.0}. 

time  propagates 

dr/dx  =  l/vh  (5.11) 

where  Vh  is  the  harmonic  mean  velocity.  For  regions  distant  from  the  source  region, 
the  mean  travel  time  propagates 

dr/dx  =  l/va  (5.12) 

where  va  is  the  arithmetic  mean  velocity,  and  equal  in  value  to  Vb,  the  bulk  filtration 
velocity.  The  harmonic  mean  velocity  is  necessarily  less  than  or  equal  to  the  arith- 
metic mean  velocity,  which  in  turn  indicates  that  the  mean  travel  time  propagates  as 
a  slower  rate  near  the  source. 

Additionally,  the  correlation  of  the  Lagrangian  velocity  persists  over  greater  dis- 
tances than  that  of  the  underlying  conductivity  field.  Solute  which  is  introduced  in 
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areas  of  high  local  flux  tend  to  move  away  from  the  source  at  relatively  high  veloci- 
ties. The  net  effect  of  what  might  be  termed  a  "preferential  flow"  of  solute  in  these 
higher  velocity  streamtubes,  is  to  induce  a  "tail"  on  a  transversely-averaged  resident 
concentration  profile,  which  increases  with  variance  of  log  conductivity. 

These  effects  are  illustrated  by  comparing  concentration  profile  created  by  averag- 
ing concentrations  along  planes  normal  to  the  mean  direction  of  flow  to  concentration 
profiles  given  by  c/c0  =  exp[-kx/v]  and  c  =  /c0  exp[-/cr[:r]],  where  r  is  given  by 

t[x]  =  —  (x  +  \ink°lk  (!  -  exp[-te/Ainfc]))  (5.13) 

"a 

Here  Ain  *  is  the  correlation  length  of  the  log  conductivity  distribution  and  b  is  a  shape 
factor  for  the  conductivity  anisotropy. 

5.4    Application:  Field  Data 

Data  from  an  elaborate  experiment  conducted  at  Columbus  Air  Force  Base,  Mis- 
sissippi were  analyzed  using  the  results  dervied  in  the  previous  section.  The  objective 
the  experiment  was  to  characterize  the  natural  attenuation  of  certain  hydrocarbons 
in  groundwater  emanating  from  subsurface  nonaqueous  phase  liquid  (NAPL)  sources 
Stauffer  et  al.  [1997].  A  NAPL  hydrocarbon  source  comprised  of  decane,  benzene, 
toluene,  ethyl  benzene,  para  xylene  (p-xylene),  and  naphthalene  was  emplaced  at  the 
site  of  the  MADE  [Boggs  et  a/.,  1992]  and  MADE-2  [Maclntyre  et  al,  1993]  experi- 
ments. This  source  released  these  poorly-soluble  organic  solutes  into  the  groundwater, 
and  the  resultant  plumes  were  monitored  at  down-gradient  multilevel  sampler  loca- 
tions. The  center  of  mass  of  the  steady-state  plume  was  estimated  using  a  trapezoidal 
rule  integration  of  concentration  data  in  space  and  by  a  simple  averaging  of  point- wise 
observations  of  concentration  via  Equation  (5.7).  Relative  degradation  and  absolute 
degradation  rates  of  the  different  constituents  are  estimated  using  Equation  (5.9). 
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5.4.1    Results  and  Discussion:  Center  of  Mass 

The  center  of  mass  of  the  steady  state  plumes  emanating  from  the  emplaced 
contaminant  source  was  estimated  by  two  methods.  The  first  might  be  considered 
a  traditional  approach,  namely  simple  trapezoidal  rule  integrations  of  the  data.  A 
known  initial  concentration  for  each  species  was  assumed  known.  This  concentration 
was  calculated  from  initial  source  mass  fractions  and  aqueous  solubilities  given  in 
Stauffer  et  al.  [1997]  and  assuming  Raoult's  law  (see  Table  5.1).  The  second  approach 
was  to  use  Equation  (5.7)  (see  Table  5.1).  The  same  initial  concentration  used  in  the 
spatial  integration  was  used  with  the  estimator  approach.  Displacements  were  taken 
as  the  distance  from  the  multilevel  sampler  to  the  source  zone. 

benzene    toluene    ethyl  benzene    p-xylene  naphthalene 
integration       4l)  21)  3l  2^3  3^8 

v/k  5.3  1.7  2.5  1.9  3.1 

Table  5.1:  The  center  of  mass  in  meters  of  steady  plumes  estimated  by  traditional 
trapezoidal  integration  of  spatial  concentration  data  and  by  the  mean  v/k  estimator 
Equation  (5.7). 

Equation  (5.7)  gives  results  similar  to  those  of  the  trapezoidal  rule  integration, 
while  offering  some  compelling  computational  advantages.  While  the  spatial  integra- 
tion is  conceptually  simple,  it  is  computationally  demanding  for  unequally  spaced 
data  in  that  the  different  "trapezoids"  must  be  calculated  and  summed  for  each 
integration.  New  data  require  recalculation  of  these  trapezoids.  Of  course,  these 
calculations  are  trivial  for  computer  codes,  but  tedious  for  "back  of  the  envelope" 
calculations.  Equation  (5.7)  is  a  fairly  simple  calculation  for  any  scientific  calculator, 
and  provides  an  excellent  tool  for  a  rapid  estimate  of  the  center  of  mass  from  sparse 
and  scattered  local  concentration.  New  data  are  easily  incorporated  into  existing  sets 
without  the  need  of  a  global  recalculation. 
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5.4.2    Results  and  Discussion:  Degradation 

The  relative  degradation  and  absolute  degradation  of  the  different  contaminants 
were  estimated  using  Equation  (5.9)  and  degradation  rates  reported  by  Maclntyre 
et  al.  [1993].  These  reported  benzene,  p-xylene,  and  naphthalene  degradation  rates 
were  estimated  from  another  experiment  conducted  at  the  site  in  which  a  "cocktail"  of 
hydrocarbons  was  injected  into  the  groundwater  and  the  fate  of  the  resultant  plume 
was  monitored  (see  Table  5.2).  This  experiment  was  conducted  before  the  NAPL 

benzene  toluene  ethyl  benzene  p-xylene 
relative        087         3A              L5  2l) 
estimated     0.0056      0.022          0.0094  0.013 
reported      0.0070       NR             NR  0.011 

Table  5.2:  Degradation  rates  relative  to  naphthalene  and  absolute  degradation  rates 
in  reciprocal  days  as  estimated  by  Equation  (5.9)  and  as  reported  by  Maclntyre  et  al. 
[1993].  NR  denotes  results  not  reported  in  Maclntyre  et  al.  [1993].  Naphthalene 
degradation  reported  to  be  0.0064 

source  was  emplaced,  and  is  in  that  sense  independent  of  the  concentration  data  as- 
sociated with  the  natural  attenuation  experiment.  Never-the-less,  the  rates  estimated 
from  the  natural  attenuation  data  were  similar  to  those  reported  by  Maclntyre  et  al. 
[1993].  The  technique  used  here  is  quite  simple,  and  does  not  require  regressions  or 
nonlinear  fits  to  transport  models. 

5.5    General  Conclusions 

These  estimators,  while  fairly  simple  and  apparently  robust,  are  limited  in  appli- 
cability to  quasi-steady  plumes,  and  it  is  not  clear  the  extent  to  which  such  plumes 
exist.  Sites  where  poorly  mobile  NAPLs  have  been  known  to  exist  for  long  times  are 
likely  candidates,  especially  if  decay  rates  are  known  to  be  large  in  comparsion  to 
transport  time  scales. 
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Local  dispersion  was  neglected,  and  local  dispersion  certainly  has  a  significant 
effect  on  the  local  concentration  variability.  However,  its  neglect  does  not  appear  to 
diminish  the  results  here  to  the  point  of  uselessness. 

The  nonstationarity  of  equal-area  Lagrangian  trajectories  in  a  steady  flow  field 
was  neglected.  For  highly  heterogeneous  conductivity  fields,  low  velocity  areas  are 
disproportionately  sampled  with  respect  to  the  flow  average  velocity.  For  a  spatially 
uniform  decay  process,  this  would  tend  to  underestimate  the  plume  extent,  since  the 
"observed  values"  will  tend  to  be  associated  with  long  travel  times.  It  should  be 
noted  that  this  sampling  bias  is  intrinsic  to  multilevel  samplers. 

The  correlation  of  velocity  and  decay  coefficient  was  neglected.  A  strong  negative 
correlation  could  result  in  a  significant  underprediction  of  the  center  of  mass  in  a 
highly  heterogeneous  field.  Again,  it  should  be  noted  that  this  is  partially  due  to  the 
sampling  bias  assumed  by  the  equal-weight  multilevel  sampler  sampling  strategy  that 
will  preferentially  sample  low  velocity  areas  with  respect  to  the  flow  weighted  mean. 


CHAPTER  6 
GENERAL  CONCLUSIONS 

Any  discipline  in  which  you  can  publish  a  paper  on  boundary  conditions 
is  an  immature  discipline.  Sten  Berglund1 

Many  of  the  results  of  this  study  are  quantitatively  and  qualitatively  the  same 
as  those  of  previous  studies.  For  instance,  chemical  engineers  have  long  known  that 
the  "mean  residence  time"  of  the  "swept  volume"  of  a  reactor  for  a  flux  injection  is 
L/U,  where  L  is  the  reactor  length  and  U  is  the  mean  filtration  velocity.  However, 
hydrologists  seek  to  understand  the  effects  of  the  heterogeneity  of  the  "reactor"  itself 
and  the  interplay  of  the  processes  as  they  manifest  themselves  "between  x  =  0  and 
x  =  L." 

An  old  premise  served  as  the  conceptual  foundation  of  this  work,  namely  that 
a  system,  in  this  case  a  heterogeneous  flow  field,  may  be  observed  in  two  ways: 
along  "straight-line,"  or  Eulerian,  trajectories  and  along  "natural"  or  Lagrangian 
trajectories  along  hydrodynamic  streamlines.  The  flow  field  was  decomposed  into  an 
aggregation  of  elements,  either  as  area  elements  that  lie  between  Eulerian  trajectories 
or  hydrodynamic  streamtubes  that  lie  in  between  streamlines.  Two  criteria  defined 
four  different  collections  of  the  flow  field  elements,  namely  that  at  some  definition 
plane,  each  element  has  either  an  equal  area  through  which  water  enters  the  element, 
or  an  equal  volume  of  water  enters  the  element  per  unit  time.   The  imaginative 


1  From  a  conversation  with  the  author  about  the  acceptance  of  Demmy  et  al.  [1999] 
for  publication. 
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reader  might  contrive  a  plethora  of  trajectory  and  trajectory  weighting  schemes.  For 
compactness,  the  following  notation  summarizes  the  four  different  trajectories 
AW/ET  equal-area  weighted  Eulerian,  or  straight  line,  trajectory 
AW/LT  equal-area  weighted  Lagrangian,  or  streamline,  trajectory 
FW/ET  equal-flow  weighted  Eulerian,  or  straight  line,  trajectory 
FW/LT  equal-flow  weighted  Lagrangian,  or  streamline,  trajectory 

Myriad  methods  exist  for  mapping  properties  observed  along  the  different  trajec- 
tories. That  is,  the  value  of  some  property  may  be  observed  along  the  trajectory  and 
recorded  at  different  intervals  in  space  or  time.  A  traditional  method  is  a  regular  Eule- 
rian gridding  that  divides  the  flow  field  into  equally-spaced  trajectories  and  "records" 
information  at  regular  intervals  along  the  trajectory.  Two  specific  parameterizations 
are  considered:  equal  intervals  in  displacement  along  the  mean  flow  direction  and 
equal  intervals  in  advective  travel  time. 

Employing  these  concepts,  the  behavior  of  different  flow  field  properties  was  exam- 
ined. These  properties,  such  as  log  conductivity  and  local  velocity,  were  "observed" 
along  the  different  trajectories  at  "control  planes,"  or  planes  equally  spaced  in  the 
mean  flow  direction,  or  "control  times,"  or  at  the  displacement  that  an  indivisible 
fluid  parcel  would  be  after  traveling  for  some  reference  time.  For  a  steady  and  irrota- 
tional  flow  field  resulting  from  a  constant  mean  gradient  applied  across  a  stationary 
Eulerian  log  conductivity  field  of  uniform  and  constant  mobile  water  content,  the 
trajectories  for  which  the  statistics  certain  properties  appeared  to  be  second-order 
stationary  were  identified.  That  is,  the  mean  and  variance  of  the  property  appeared 
to  have  a  constant  value  at  all  observation  points  along  the  trajectory.  The  log  con- 
ductivity and  velocity  observed  at  control  planes  appeared  stationary  for  the  equally- 
spaced  Eulerian  trajectories  (AW/ET)  and  Lagrangian  trajectories  separated  by  an 
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equal  discharge  (FW/LT).  For  equally-spaced  control  times,  only  Lagrangian  trajec- 
tories equally-spaced  at  a  reference  plane  (AW/LT)  appeared  to  exhibit  stationary 
log  conductivity  and  velocity  statistics. 

The  stationarity  of  the  properties  greatly  simplifies  the  development  of  predictors 
of  solute  transport.  When  taking  the  expected  value  of  the  integral  of  some  property 
along  a  trajectory,  say  the  integral  of  velocity  in  time,  the  time-varying  velocity 
function  may  be  replaced  with  its  time  average. 

As  implied  by  the  stationarity  of  the  statistics  along  some  trajectory  collections, 
but  not  others,  the  numerical  statistics  of  properties  varied  between  the  different 
collections.  Adopting  simple  and  accurate  models  for  the  statistical  properties  of 
some  property  in  a  reference  collection,  the  statistics  of  the  properties  along  other 
collections,  including  the  observed  nonstationarities,  were  well  predicted.  This  had 
significant  application  in  the  evaluation  of  the  effect  of  injection  mode  upon  solute 
transport. 

Two  methods  by  which  solute  might  be  uniformly  introduced  into  a  system  were 
considered.  The  first  was  a  resident  injection  were  a  fixed  volume  at  the  system  inlet 
was  uniformly  filled  with  solute  at  some  reference  concentration.  The  second  was  to 
allow  influent  water  at  some  reference  concentration  to  carry  the  mass  into  the  system. 
The  displacement  and  travel  time  statistics  of  AW/LT  correspond  to  solute  transport 
associated  with  a  uniform  resident  injection.  The  transport  of  solute  injected  in  flux 
is  described  by  the  displacement  and  travel  time  statistics  of  the  FW/LT. 

The  so-called  consistent  first-order  approximations  of  the  inverse  velocity  were 
shown  to  be  built  upon  a  flawed  premise,  but  are  none-the-less  quantitatively  cor- 
rect and  robust.  The  premise  is  that  the  Lagrangian  flow  field  may  be  accurately 
approximated  by  the  Eulerian  flow  field.  While  this  may  be  true  for  the  weakest 
of  heterogeneities,  this  approximation  breaks  down  rapidly  with  increasing  hetero- 
geneity if  it  is  consistently  applied.  As  inconsistently  formulated,  the  approximation 
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<  1/v  >~  1/U  is  quantitatively  equal  to  the  inverse  of  the  harmonic  mean  of  the 
Lagrangian  velocity  <  1/v  >=  l/Vh.  Moreover,  the  time-average  Lagrangian  veloc- 
ity is  quantitatively  equal  to  the  space  average  Eulerian  velocity,  and  these  equalities 
contribute  to  the  robustness  of  other  approximations.  Were  the  Lagrangian  field  to 
be  consistently  approximated  by  the  Eulerian,  travel  times  would  be  greatly  over- 
predicted  and  displacements  greatly  under-predicted.  It  is,  in  fact,  the  consistent 
inconsistent  approximation  of  the  Lagrangian  by  the  Eulerian  that  contributes  to 
much  of  the  observed  robustness  of  the  Lagrangian  transport  theory.  While  not 
many  "new"  equations  were  proposed  or  so  many  new  phenomena  were  unearthed, 
it  is  hoped  that  a  modicum  of  understanding  was  fostered.  This  is  the  primary 
contribution  of  this  work. 

These  concepts  were  illustrated  with  an  analysis  of  the  transport  of  an  solute  ex- 
periencing a  linear  equilibrium  sorption  in  a  field  characterized  by  a  variable  sorption 
coefficient.  While  the  process  has  been  studied  before,  it  is  important  to  illustrate 
the  new  concepts  and  approaches  here  with  familiar  examples.  The  sorption  study 
was  restricted  to  transport  as  characterized  by  mass  breakthrough  at  control  planes. 
An  interesting  new  result  is  that  injection  mode  strongly  affects  the  propagation  of 
the  mean  breakthrough  time.  A  uniform  resident  injection  exhibits  a  nonlinearity 
in  the  mean,  and  this  nonlinearity  is  predictably  enhanced  by  negative  correlation 
between  the  sorption  coefficient  and  the  underlying  log  conductivity  field.  Injection 
in  flux,  however,  results  in  a  linear  propagation  of  mean  arrival  time  in  space,  and 
mitigates  the  effect  upon  correlation.  Arrival  time  variance,  however,  is  dominated  by 
the  correlation  between  sorption  coefficient  and  the  log  conductivity  field.  Injection 
mode  has  comparatively  little  effect.  This  has  profound  implications  for  tracer  test 
analyses  based  upon  the  arrival  time  variance  of  breakthrough  curves. 

The  concepts  developed  in  this  work  are  further  illustrated  with  the  develop- 
ment and  analysis  of  an  estimator  of  the  center  of  mass  of  a  plume  resulting  from 
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continuously  injected  solute  subject  to  a  first-order  decay.  The  assumption  of  non- 
interacting  streamtubes  with  constant  "effective"  properties  allows  solution  of  the 
advection-decay  equation  for  a  single  streamtube.  An  ergodic-like  hypothesis  al- 
lows this  solution  to  a  distribution  of  streamtubes.  This  is  what  might  be  termed 
a  "classic"  streamtube  approach.  While  surprisingly  robust  and  perhaps  worthy  of 
application  in  field  settings,  the  estimator  shows  a  bias  explained  using  the  trajectory 
based  analysis  developed  in  this  work. 

One  measure  of  quality  of  a  work  is  the  number  of  questions  that  it  answers.  Per- 
haps another  is  the  number  of  interesting  questions  that  it  begs.  Several  issues  are 
left  for  further  study.  A  few  specific  questions  are  addressed  first,  then  a  few  more 
general.  The  relationship  between  an  Eulerian  "reference"  system  and  its  Lagrangian 
counterpart  remains  incompletely  understood.  A  primary  area  that  requires  attention 
is  the  increasing  Lagrangian  velocity  correlation  length  with  increasing  log  conduc- 
tivity variance.  A  great  deal  of  work  remains  for  predicting  the  spatial  moments  of  a 
reactive  solute. 

Numerical  constraints  limited  this  study  to  relatively  mild  heterogeneity.  How- 
ever, a  highly  heterogeneous  medium,  might  be  replaceable  by  two  or  more  media 
characterized  by  lesser  heterogeneity,  and  a  systematic  "replacement  system"  would 
greatly  enhance  computational  efforts.  The  turning  bands  method  employed  in  the 
generation  of  the  random  fields  was  computationally  expensive,  in  comparison  to, 
say,  the  particle  tracking  algorithms  or  data  analysis  programs,  and  to  the  flow  solver 
for  weak  heterogeneity.  The  method  requires  the  generation  of  several  independent 
normal  line  processes.  This  problem  is  well-suited  to  parallel  and  distributed  compu- 
tational techniques,  and  an  easy-to-use  implementation  would  certainly  find  wide  use. 
Several  research  groups  are  working  on  parallel  and  distributed  flow  codes,  and  the 
research  community  would  benefit  greatly  from  a  wider  dissemination  of  or  greater 
access  to  these  codes. 
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Thinking  more  broadly,  some  of  the  concepts  employed  here  might  find  application 
in  other  fields.  The  flow  of  water  through  a  wetland  system  might  be  reasonably 
approximated  as  a  two-dimensional  aquifer-like  object.  Transport  of  contaminants 
might  be  modeled  using  similar  techniques,  especially  if  the  transport  time  scales 
are  relatively  short  compared  to  system  variability.  In  general,  heterogeneity  will 
imply  some  sort  of  "preferential  flow,"  and  the  identification  of  preferential  flows  in 
heterogeneous  systems  should  be  paramount  consideration  in  the  design  of  monitoring 
systems. 


APPENDIX  A 
EULERIAN  VELOCITY  DISTRIBUTION  PROPERTIES 

Consider  lognormally  distributed  random  variable  V  with  moments  <  V  >=  U 

and  var(V)  =  S.  These  moments  are  related  to  the  parameters  \i  and  a2  via 

U  =  exp(At  +  <t2/2)  (A.l) 

S  =  exp(2/x  +  2a2)  -  exp(2//  +  a2)  (A.2) 
Let  M  =  exp(/i)  and  Z  =  exp(cr2).  Rewrite  the  previous  equations  as 

U  =  MZ1/2  (A.3) 


S  =  M2(Z2  -  Z)  (A.4) 


Solving  for  M  and  Z  yields 


U2 

M  -  TfTF  (A'5) 

(A.6) 

Expanding  the  definitions  of  M  and  Z  and  taking  the  natural  log  of  both  sides  of  the 
previous  two  equations  yields 

"2  =  In  (2±£)  (A.8) 


125 


126 

The  effective  conductivity  of  a  porous  medium  may  be  denned  as  the  numerical 
conductivity  value  that  satisfies  the  Darcy  relationship 

Q  =  KeJ  (A.9) 

where  Q  is  the  specific  discharge  and  J  is  the  difference  in  hydraulic  head  across  the 
formation.  There  are  several  expressions  for  the  effective  conductivity  of  a  stationary 
heterogenous  conductivity  fields  (see  Gelhar  [1993]).  The  general  form  of  the  effective 
conductivity  is 

Ke  =  Kgf(alk,g)  (A.10) 

where  Kg  and  afnk  are  the  geometric  mean  and  the  variance  of  the  log  conductivity 
field,  respectively,  g  is  a  function  of  aquifer  anisotropy,  and  /  is  a  function  that  scales 
the  geometric  mean  to  the  proper  effective  value.  For  a  two-dimensional  isotropic 
aquifer,  f{cr2nk)  =  1.  Thus,  for  an  aquifer  of  homogeneous  mobile  water  content  6, 
an  estimate  of  the  spatial  mean  seepage  velocity  U  is  given  by 

U=IY  (A.11) 

The  point  velocity  variance  for  such  a  field  is  S  =  bU2a2nk,  where  b  =  3/8. 
Substitution  of  Equation  (A. 11)  into  this  expression  yields 

bK2J2o2 

S=      9 Q2  (A.12) 


Substitution  of  Equations  (A. 11)  and  (A.12)  into  Equations  (A. 7)  and  (A. 8) 

0y/TTbalk 

a2  =  ln(l  +  ba2nk)  (AAA) 


I*  =  ]n  \  —j£sL=  ]  (A.13) 


APPENDIX  B 
VELOCITY  -  LOG  CONDUCTIVITY  CORRELATION 


A  simplified  correlation  relationship  between  the  velocity  field  and  the  log  con- 
ductivity is  sought.  Graham  and  McLaughlin  [1989b]  and  Rubin  and  Dagan  [1992] 
derived  expressions  based  upon  perturbation  expansions  of  Darcy's  law  and  spectral 
methods.  Graham  and  McLaughlin  [1989b]  assume  a  "hole-type"  covariance  function 
for  the  log  conductivity  field,  in  order  to  assure  a  complementary  stationary  head  gra- 
dient field.  Rubin  and  Dagan  [1992]  assume  an  exponentially  correlated  anisotropic 
log  conductivity  field. 

The  velocity-log  conductivity  covariance  functions  for  both  analytical  cases  shows 
qualitatively  similar  behavior  (see  Figure  B.l),  and  an  simplified  estimator  would  be 
applicable  for  both  log  conductivity  correlation  models.  Although  these  rigorously- 
developed  estimators  are  attractive  theoretically,  the  mundane  and  oft-neglected  task 
of  mapping  theoretical  results  into  application  often  requires  simplifying  assumptions. 
Presented  for  your  consideration  is  the  covariance  function  of  Graham  and  McLaughlin 
[1989b]  (their  Equation  (B5) 


where  a  =  7r/(4A/)  and  K{  is  the  modified  Bessel  function  of  order  i.  The  integral 
scale  of  this  covariance  function  is  A/.  An  exponential  approximation  to  this  function 
that  reproduces  the  integral  scale,  general  trend,  and  point  covariance  of  that  given 
by  Graham  and  McLaughlin  [1989b]  is  sought.  A  function  of  the  form 


-(a/2)  (C#i[aC]-aC2ffo[aC])} 


(B.l) 


29 


exp[-x/(2A)] 


(B.2) 
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Figure  B.l:  Comparison  of  velocity  and  log  conductivity  correlation  functions  from 
Graham  and  McLaughlin  [1989b]  and  Rubin  and  Dagan  [1992]  and  a  simplified  expo- 
nential model.  All  functions  are  normalized  by  KgJafnk/9.  The  apparent  "roughness" 
of  the  Rubin  and  Dagan  [1992]  function  is  due  to  estimation  errors  associated  with 
inferring  their  function  from  graphical  data. 
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meets  these  criteria  (see  Figure  B.l). 

The  covariance  function  presented  by  Rubin  and  Dagan  [1992]  is  the  sum  of 
three  two-function  products  that  involve  integrations  of  Bessel  functions.  For  two- 
dimensional  porous  systems,  the  simplified  function  condenses  that  monumental  work 
into  an  elegant  and  concise  form. 

The  simplified  expression  is  tested  against  correlation  functions  estimated  from 
the  velocity  and  log  conductivity  data  taken  from  the  simulations  (see  Figure  B.2). 
The  simplified  estimator  gives  a  reasonable  a  priori  estimate  of  the  velocity  and 
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Figure  B.2:  Comparison  of  the  simplified  velocity-log  conductivity  correlation  func- 
tion to  simulation  data  for  aj  =  {1/2, 1}. 

log  conductivity  correlation  observed  along  the  equal-flux  Lagrangian  trajectories. 
The  reasonable  match  to  analytical  expressions  and  the  observed  data  enhance  our 
confidence  in  this  simplified  estimator. 

The  results  of  Rubin  and  Dagan  [1992]  (their  Figure  7)  indicate  that  anisotropy 
decreases  the  point  covariance  and  increases  the  long  range  correlation  effects.  Yet 
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the  functional  form  appears  to  remain  exponential.  Therefore,  it  is  suggested  that 
future  work  focus  upon  generalizing  the  simple  results  found  here  to  three-dimensions 
and  general  anisotropy. 


APPENDIX  C 
CONSISTENT  FIRST  ORDER  APPROXIMATION 


C.l    Travel  Time 


Usage  of  a  consistent  first-order  approximation  for  travel  time  permeates  the 
stochastic  Lagrangian  literature.  This  approximation  is  quantitatively  correct  for 
collections  of  equal-flux  streamtubes.  However,  the  premise  upon  which  it  is  built  is 
incorrect.  Consider  a  stationary  lognormal  process  V.  The  lognormal  distribution  is 
selected  for  convenience,  as  a  random  variable  with  values  guaranteed  to  be  greater 
than  zero  is  required.  Since  V  is  stationary  by  definition,  then  so  is  1/V. 

Consider  the  expected  value  of  the  integral 


The  order  of  evaluation  may  be  interchanged,  taking  the  expected  value  of  the  argu- 
ment of  the  integral  thus: 


where  Vh  is  the  point  harmonic  mean  of  the  random  process  V.  Consider  now  a 
perturbation  expression  for  the  random  process  V  =  V(l  +  v)  where  V  is  the  mean  of 
V,  and  v  is  a  zero-mean  random  process  representing  the  perturbation  of  the  process 
about  its  mean.  Inserting  this  expression  into  Equation  (C.l)  yields 


(C.l) 


DO 


<dx/V[x]> 


(C.2) 


=  x/Vh 


(C.3) 
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Expand  (1  +  v)  in  a  series, 

1  f°° 

<r\x)>=Tj     <{1-v[x]  +  v[x]2  )>dx  (C.4) 

v  Jo 

Should  this  expression  be  truncated  at  first  order,  the  leading  term  <  r  x/V  is  left. 
However,  from  Equation  (C.2),  the  exact  expression  is  <  r  >=  x/Vh.  <  r  >=  x/V 
is  a  well-established  result  in  the  literature  for  the  first  temporal  moment  of  a  solute 
plume  injected  in  flux  Kreft  and  Zuber  [1978]. 

The  answer  to  this  riddle  lies  in  the  examination  of  the  properties  of  the  flow 
field.  Even  if  the  deviations  are  relatively  small  from  the  straight-line,  or  Eulerian, 
trajectory,  the  statistics  of  the  properties  recorded  along  the  actual,  or  Lagrangian, 
trajectory  can  be  quite  different.  For  instance,  the  mean  velocity  recorded  along 
equal-flow  weight  streamlines  in  an  quasi-infinite  constant  mean  gradient  flow  in 
a  exponentially  correlated  lognormal  log  conductivity  field  is  approximately  V  = 
(1  +  bafnk)U  where  U  is  the  spatially  averaged  Eulerian  velocity.  The  variance  of 
the  the  velocity  perturbation  along  these  trajectories  is  approximately  o\  =  bafnk. 
Evaluating  the  expected  value  in  Equation  (C.4)  and  retaining  the  second  order  term 
yields 


<r[x)>=t(l+ol)  (C.5) 
Substituting  our  approximate  values  into  the  preceding  equation  yields 

\l  +  h(JLk)u  (C.6) 
=  x/U 

The  so-called  small-deviation  assumption  is  a  poor  one.  The  approximation  of  the 
Lagrangian  field  by  the  Eulerian  field  works  because  the  correlation  properties  are 
quite  similar  and  the  "coincidence"  that  the  harmonic  mean  Lagrangian  velocity  is 
approximately  equal  to  the  arithmetic  mean  Eulerian  velocity. 
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C.2  Displacements 

The  center  of  mass  of  a  solute  plume  uniformly  injected  in  space  is  known  to 
move  <x[f]>=  Ut,  where  U,  again,  is  the  spatially  averaged  Eulerian  velocity.  If  the 
solute  parcels  move  along  Lagrangian  trajectories,  and  the  statistics  of  Lagrangian 
and  Eulerian  fields  are  known  to  be  quite  different,  why  does  this  work? 

The  answer  is  that  the  temporally  averaged  velocity  of  equal-area  Lagrangian 
streamtubes  is  stationary  in  time  and  equal  to  the  spatially  averaged  Eulerian  ve- 
locity. In  fact,  the  temporally  averaged  Eulerian  velocity  is  lower  than  its  spatially 
averaged  counterpart.  So  once  again,  the  small-displacement  assumption,  if  consis- 
tently applied,  would  result  in  a  very  poor  result.  What  is  remarkable  is  that  the 
correlation  properties  of  the  two  fields  are  so  similar,  and  herein  lies  the  strength  of 
the  approximations  of  Lagrangian  fields  with  Eulerian  fields. 
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